{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Random Signals and LTI-Systems\n",
    "\n",
    "*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Auto-Correlation Function\n",
    "\n",
    "The auto-correlation function (ACF) $\\varphi_{yy}[\\kappa]$ of the output signal of an LTI system $y[k] = \\mathcal{H} \\{ x[k] \\}$ is derived. It is assumed that the input signal is a wide-sense stationary (WSS) real-valued random process and that the LTI system has a real-valued impulse repsonse $h[k] \\in \\mathbb{R}$. \n",
    "\n",
    "Introducing the output relation $y[k] = h[k] * x[k]$ of an LTI system into the definition of the ACF and rearranging terms yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\varphi_{yy}[\\kappa] &= E \\{ y[k+\\kappa] \\cdot y[k] \\} \\\\\n",
    "&= E \\left\\{  \\sum_{\\mu = -\\infty}^{\\infty} h[\\mu] \\; x[k+\\kappa-\\mu] \\cdot \n",
    "\\sum_{\\nu = -\\infty}^{\\infty} h[\\nu] \\; x[k-\\nu] \\right\\} \\\\\n",
    "&= \\underbrace{h[\\kappa] * h[-\\kappa]}_{\\varphi_{hh}[\\kappa]} * \\varphi_{xx}[\\kappa]\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "where the ACF $\\varphi_{hh}[\\kappa]$ of the deterministic impulse response $h[k]$ is commonly termed as *filter ACF*. This is related to the [link between ACF and convolution](../random_signals/correlation_functions.ipynb#Definition). The relation above is known as the *Wiener-Lee theorem*. It states that the ACF of the output $\\varphi_{yy}[\\kappa]$ of an LTI system is given by the convolution of the input signal's ACF $\\varphi_{xx}[\\kappa]$ with the filter ACF $\\varphi_{hh}[\\kappa]$. For a system which just attenuates the input signal $y[k] = A \\cdot x[k]$ with $A \\in \\mathbb{R}$, the ACF at the output is given as $\\varphi_{yy}[\\kappa] = A^2 \\cdot \\varphi_{xx}[\\kappa]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - System Response to White Noise\n",
    "\n",
    "Let's assume that the wide-sense ergodic input signal $x[k]$ of an LTI system with impulse response $h[k] = \\text{rect}_N[k]$ is normal distributed white noise. Introducing $\\varphi_{xx}[\\kappa] = N_0\\, \\delta[\\kappa]$ and $h[k]$ into the Wiener-Lee theorem yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{yy}[\\kappa] = N_0 \\cdot \\varphi_{hh}[\\kappa] = N_0 \\cdot (\\text{rect}_N[\\kappa] * \\text{rect}_N[-\\kappa])\n",
    "\\end{equation}\n",
    "\n",
    "The example is evaluated numerically for $N_0 = 1$ and $N=5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAGFCAYAAABT4e8GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZRddX3v8feXSYAR0GkkIhlAUXGEGjRWKSx764DaYH0gom19qNrWXvRWWr21UaK1re1VaLPa67I+VFofelsfQA3R1mqkxVlWFw+KQYLiVAQKTlBEHCAwQDL53j/OHjhMZubskzlzfjNz3q+1ZiVz9u/s/Z3f2Xufz9n7t/eJzESSJEnlHFC6AEmSpF5nIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJi1xEfCcihkvXMZOI+FhE/J/SdbQSEUMRsT0i7oqIPyhdz0Lo1noSETdGxHP243kZEXdHxLumPX5zRKybof0lEXFvRHxtPvVKS4WBTFog1RvXRETsavp5X43nPOTNLjN/PjNHFqi+tt9Y92M5IxHxs4g4aIZpr4iIb1Z9c0tEfDEifmlajdP7cM1+lPEWYCQzD8vM987n75muk/04n3kt1HrSYU/JzLdP/RIRPwccCVw7vWFmnga8vou1SUUZyKSF9cLMPLTp5+zSBXVTRDwW+B9AAi+aNu0PgfcA7waOAI4BPgCcMW020/tw536U8hjgO/vxPC2stcB1mXlv6UKk0gxkUgER8daIGKtOoY1GxLMj4p9ohJJ/qY4EvaVq+8BRk+r/GyPi6ur0z4cj4ojqyNJdEfHv1VGHqeWcExE/qKZ9NyJeXD0+27LWRMRnI+InEXHD9NN7EbEuIr5Vze8C4OAWf+qrgcuAjwGvaZrPI4A/B96QmVsy8+7M3J2Z/5KZG/ejP4+vjsSNV6fuXtQ07RLgVOB91d/6xDafnxHxhKbfHzhNO1M/Vq/Rpqq/fxYRH42Ig1vNb7bXZIZa91l3qscfcnQtIp7WdJr20xFxQVPdN0bEH1Xr0R3VtOYaZ1xvWrwGh0bEZEQc2fTYk6sjn4fN8rQTgWuqtg+LiE9ExJaIOLTV8qTlxkAmdVlEDAFnA8/IzMOA9cCNmfkq4CYePCL0V7PM4iXAc4EnAi8Evgi8DTicxjbdHKJ+QOMI1SOAdwL/HBFHzrSsiDgA+Bfg28Ag8GzgTRGxvqr7QGAr8E/AKuDTVS1zeTXw8epnfUQcUT1+Co0wd1GL57cUESurur8MPAr4feDjVT9Pnfr6T+Ds6m/9r3aeP5c5XrNX0nhdH0/jdfrjecyrudYZ150Z2h1Io28/RuO1+iQwPVT9OnA6cCyNYPRbTdNmXG9a1L8L+B7wtKaHzwPenZl3zfK0E4EdEXEs8DVgFHhJNS+ppxjIpIW1tTrqMvXzP4FJ4CDghIhYmZk3ZuYP2pjn32bmjzNzjEbQuDwzt2fmfTTehB8YIJ2Zn87MnZm5NzMvAL4PnDTLfJ8BrM7MP8/M+zPzeuDvgZdV008GVgLvqY5mfQb4xmxFRmMs2GOACzPzShpv8q+oJj8SuC0z99T4e5v7cOsM008GDgXOq+q+BPhX4OU15t2J58/kfZl5c2beDrxrnvNqVnfdORlYAby3eq22AFdMa/Peat24nUYgferUhDbXm2bfoApkEfHLwAnAh+Zov5bGGLJLgHdm5jszM2ssR1p2DGTSwtqQmQNNP3+fmdcBbwL+DLg1Ij4V7Q1U/3HT/ydm+P2B0z0R8eqIuGoq0ABPpnEkbSaPAdY0B0gaR96mjmqtAcamvWH+9xx1vgb4cmbeVv3+CR48bflT4PCIWDHH86c09+GGGaavAW7OzL3T6hqsMe9OPH8mN0+b1/5ciLCPNtadmV6rm6e1+VHT/+9h/9ebZg8EMuCvgHdk5v0zNYyIqOb7YuDvMvNzNeYvLVsGMqmAzPxEZk4dQUrgL6cmdWoZEfEYGke4zgYemZkDNMbrxCzLuhm4YVqAPCwzf7WafgswWL2RTjlmlmX30zgl9qyI+FFE/Aj438BTIuIpwKXAvcBMAatdO4Gjq1OuzXWNdej59wAPa5r26GnPn+k1O3ravJovRJhrfi1f/znWnWYzvVZHz9BuHzXWm7l8A3haRLwE6KdxqnQ2x1b/Pgd4c0Q8vU590nJlIJO6LBr3xDotGreBuJfGUa3JavKPgcd1aFGH0HjD/km13N+mcURiyvRlXQHcWQ0a74+IvmpQ9jOq6ZcCe4A/iIgVEXEms5/G2kDjbzqBxqmwpwLH0zjF+urMvAP4E+D9EbGhGtC9MiKeFxGzjZ2bzeXA3cBbqnkM0xhb96kOPf8q4BVVf5wOPGva82d6zd4QEUdFxCoaRxkvaJo21/zmfP1brDvNLq0eP7t6rc6g3ilHaL3ezOXbNALmXwPnTDvqON2JwNWZuQM4C7io1Tg1aTkzkEkLa+qKuamfi2iMAToPuI3GaaNH0XjTBjgX+OPqVNEfzWfBmfldGm+Ml9J4o18LfL2pyUOWlZmTNILIU4Ebqvr+gcbAbqpTT2fSGPz9M+A3gC2zLP41wEcz86bM/NHUD/A+4JURsSIz/wb4QxoD3n9C4wjd2TQuHGjn77yfxi01nlfV/AEaoe97HXr+G2n0yziNwfrT65vpNfsEjYsErq9+mm+eO9f8Wr3+c6070/+mM4HXVsv5TRrj4u6bqy+q57Zab+Z67n3ADhoXqXyxRfO1wNXV87YC59MYL9jqyl1pWQrHT0pS50TEjcDvZua/l66lWURcTmOs1kcXcBkHAtcBv56Zl02bdi+NQPjezHxHjXldTOPihCsy89kLUa+0mNQZUCtJWmIi4lk0biNxG40jcScCX1rgxf4p8PXpYQwgM9s68pWZz+1YVdISYCCTpOVpCLiQxtWTPwBempm3LMSCIuJpwFdonIJseRNZSfvylKUkSVJhDuqXJEkqzEAmSZJU2JIYQzYwMJBPeMITWjdUx9x9990ccsghpcvoKfZ599nn3Wefd5993n1XXnnlbZm5up3nLIlAdsQRR/DNb36zdBk9ZWRkhOHh4dJl9BT7vPvs8+6zz7vPPu++iJjra+Vm5ClLSZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpsBWlFhwRNwJ3AZPAnsx8eqlaJEmSSioWyCqnZuZthWuQ1GO2bh9j87ZRxsYnGLzsEjauH2LDusHSZUnqYaUDmSR11dbtY2zasoOJ3ZMAjI1PsGnLDgBDmaRiIjPLLDjiBuBnQAIfyszzp00/CzgLYPXq1b9w4YUXdr/IHrZr1y4OPfTQ0mX0FPu8O948cg8/vXff/d4jDw7+evhhBSrqLa7n3Wefd9+pp556ZbtDsUoGsjWZuTMiHgVcDPx+Zn51prZDQ0M5Ojra3QJ73MjICMPDw6XL6Cn2eXcce84XmGmvF8AN5z2/2+X0HNfz7rPPuy8i2g5kxa6yzMyd1b+3AhcBJ5WqRVLvWDPQ39bjktQNRQJZRBwSEYdN/R/4FeCaErVI6i0b1w/Rv7LvIY/1r+xj4/qhQhVJUrlB/UcAF0XEVA2fyMwvFapFUg+ZGrj/ls9czf2Texkc6PcqS0nFFQlkmXk98JQSy5akDesG+eQVNzE+Ps62t55WuhxJ8k79kiRJpRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJha0oteCI6AO+CYxl5gtK1SFpedm6fYzN20bZOT7BmoF+Nq4fYsO6weLzkqS5FAtkwBuBa4GHF6xB0jKydfsYm7bsYGL3JABj4xNs2rIDoO0g1cl5SVIrRU5ZRsRRwPOBfyixfEnL0+Ztow8EqCkTuyfZvG206LwkqZVSY8jeA7wF2Fto+ZKWoZ3jE2093q15SVIrXT9lGREvAG7NzCsjYniOdmcBZwGsXr2akZGR7hQoAHbt2mWfd5l9Pn+rDg5+em/O+Pj0vh0fn2BycnLWPm9nXqrP9bz77POlITL33eEs6AIjzgVeBewBDqYxhmxLZv7mbM8ZGhrK0VFPE3TTyMgIw8PDpcvoKfb5/E0f9wXQv7KPc89cu8+4r9/40KWMj4+z7a3Pm/e8VJ/reffZ590XEVdm5tPbeU7XT1lm5qbMPCozHwu8DLhkrjAmSXVtWDfIuWeu5cC+xq5tcKB/vwNUJ+clSa2UvMpSkjpuw7pBPnnFTQBc8LpTFs28JGkuRQNZZo4AIyVrkCRJKs079UuSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBW2osRCI+Jg4KvAQVUNn8nMPy1Ri6SlY+v2MTZvG2Xn+ARrBvrZuH6IDesGrUnSklckkAH3Aadl5q6IWAl8LSK+mJmXFapH0iK3dfsYm7bsYGL3JABj4xNs2rIDoFgAWow1SVqaipyyzIZd1a8rq58sUYukpWHzttEHgs+Uid2TbN42WqiixVmTpKWp2BiyiOiLiKuAW4GLM/PyUrVIWvx2jk+09Xg3LMaaJC1NpU5ZkpmTwFMjYgC4KCKenJnXTE2PiLOAswBWr17NyMhImUJ71K5du+zzLrPP57bq4OCn9+57IH3VwbFPv41XgWiu/hwfn2BycrJln881r3ZqUoPreffZ50tDsUA2JTPHI2IEOB24punx84HzAYaGhnJ4eLhIfb1qZGQE+7y77PO5veMRDx2vBdC/so93nLGW4WnjtT44eikAw8OnzDq/D45eyvj4eMs+n2te7dSkBtfz7rPPl4YipywjYnV1ZIyI6AeeA3yvRC2SloYN6wY598y1HNjX2G0NDvRz7plriw6eX4w1SVqaSh0hOxL4x4jooxEKL8zMfy1Ui6QlYsO6QT55xU0AXPC62Y9+ddNirEnS0lM7kEXEqhrN9mbmeKtGmXk1sK7usiVJkpazdo6Q7ax+Yo42fcAx86pIkiSpx7QTyK7NzDmPakXE9nnWI0mS1HPaGdRfZ3CEAygkSZLaVDuQZea9ABFxxvRpEXFAcxtJkiTVtz+3vXhdRPwiPHC3/d/BW1ZIkiTtt/257cXLgc9FxBeA/wXsAF7d0aokSZJ6yP4EsmcCbwc+DvxWZo50tCJJkqQesz+B7KXAWuBw4EMRcTFwTWb+XUcrkyRJ6hFtB7LM/B2AiAjgOBrhbG2H65IkSeoZ+/3VSZmZwH9VP5/tWEWSJEk9pvZVlhHxrU60kSRJ0kO1c4Ts+Ii4eo7pATxinvVIkiT1nHYC2ZNqtJnc30IkSZJ6Ve1Alpn/PdPjEXFQZt7XuZIkSZJ6y/7cqX+6D0TEczswH0mSpJ4070CWma8FjomI90bE4R2oSZIkqafMO5BFxHrgWODxwD9ExIvnXZUkSVIP6cQpy8cCH87M52fmBuC0DsxTkiSpZ3QikJ0EPKHp97d3YJ6SJEk9o+NjyDLzzg7UJUmS1DMcQyZJklTYfn2XZUQ8AzgQ+A7waOAjmXl9Ne1vgYs6VqEkSdIy13Ygi4itwFHAzcAJwNeBTzc1cQyZJElSG1qesoyIEyLin5seejLwBuA3MnMI+BrwwamJjiGTJElqT50jZP8BnNL0+3nAW4GnRMQ9wA7gWdXd+q/KzJ90vkxJkqTlq86g/l8B3tX0+wuBCzPz8cAvAX8H9AEvA77U8QolSZKWuZZHyDJzB/DKpod+F/hoRGwCrgGGgC9Vt7+QJElSm9oe1F+dknxBRKwB1gJ3ZOZlHa9MkiSpR+zXbS8AMnMnsLODtUiSJPWkTnx1kiRJkubBQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVNiK0gVI0tbtY2zeNsrO8QnWDPSzcf0QG9YNli6rY5b73ydp/gxkkoraun2MTVt2MLF7EoCx8Qk2bdkBsCxCy3L/+yR1RpFTlhFxdER8JSKujYjvRMQbS9QhqbzN20YfCCtTJnZPsnnbaKGKOmu5/32SOqPUEbI9wJsz81sRcRhwZURcnJnfLVSPpEJ2jk+09fhSs9z/PkmdUeQIWWbekpnfqv5/F3At4LF7qQetGehv6/GlZrn/fZI6o/hVlhHxWGAdcHnZSiSVsHH9EP0r+x7yWP/KPjauHypUUWct979PUmcUHdQfEYcCnwXelJl3Tpt2FnAWwOrVqxkZGel+gT1s165d9nmX9WqfDwCvOr6Pj1wzyZ698MiDg5c8sY+BO77PyMj392k/Xp3qm6uv6raZnJxs2efzXV67f99y16vreUn2+dJQLJBFxEoaYezjmbll+vTMPB84H2BoaCiHh4e7W2CPGxkZwT7vrl7u82Hg2x+6FIALXnfKnG0/ONpoNzw8e7u6bcbHx1v2eSeWN0z9v2+56+X1vBT7fGkodZVlAB8Grs3MvylRgyRJ0mJRagzZM4FXAadFxFXVz68WqkWSJKmoIqcsM/NrQJRYtiRJ0mJT/CpLSZKkXmcgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJha0oXYCk5W3r9jE2bxtl5/gEawb62bh+iA3rBkuXtejYT1JvM5BJWjBbt4+xacsOJnZPAjA2PsGmLTsADBtN7CdJRU5ZRsRHIuLWiLimxPIldcfmbaMPhIwpE7sn2bxttFBFi5P9JKnUGLKPAacXWrakLtk5PtHW473KfpJUJJBl5leB20ssW1L3rBnob+vxXmU/SfIqS0kLZuP6IfpX9j3ksf6VfWxcP1SoosXJfpK0aAf1R8RZwFkAq1evZmRkpGxBPWbXrl32eZctxz4fAF51fB8fuWaSPXvhkQcHL3liHwN3fJ+Rke8/pO14dXquVR/UaVe3zeTkZFeXN1ubdvppqVuO6/liZ58vDYs2kGXm+cD5AENDQzk8PFy2oB4zMjKCfd5dy7XPh4Fvf+hSAC543SmztvvgaKPN8PDsbeq2q9tmfHy8ZZ93cnlztRmmXj8tdct1PV/M7POlwVOWkiRJhZW67cUngUuBoYj4YUS8tkQdkiRJi0GRU5aZ+fISy5UkSVqMPGUpSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmErSheg7tm6fYzN20bZOT7BmoF+Nq4fYsO6wRnbjI1PMHjZJTO2kSQtfnX353XeG7TwDGQ9Yuv2MTZt2cHE7kkAxsYn2LRlB8ADG16dNpKkxa/u/tz9/uLhKcsesXnb6AMb3JSJ3ZNs3jbaVpspW7eP8czzLuHYc77AM8+7hK3bxxamcEnSPlrtg+vuz9vZ72theYRsmWh1yHnn+MSMz2t+vE6bqWX5iUqSyqizD667P29nv+9pzYXlEbJlYGrjHBufIHlw42z+xLRmoH/G5zY/XqcN1P9E5VE0SWpfJ45+1d2f12lX5z1G82cgWwbqbJwb1w/Rv7LvIW36V/axcf1QW22g3icqN2BJal+dfWedfXDd/Xmddp7W7A4D2SJX5yhTnY1zw7pBzj1zLQf2NV7ywYF+zj1z7UMOOddpA/U+UbkBS1L7OnX0q+7+vE67dk5relZk/xnIFrG6R5nqHpresG6QdccM8IvHruLr55w24/n/qTZDP3fArG3qfKKquwFLkh7UyaNfdfbnze1me2/wtGZ3GMgWsbpHmepunJ1S5xNV3ZDoJypJvaLO/q6TR786xdOa3eFVlgV14spIePCqmrd85mrun9zLYBeugNmwbpBPXnETABe87pR9pm9cP/SQq4Bg3w3YqzUl9Yq6+7s6+86p58y1D+6kOu8xXq05fwayQupsnGsG+hmbYSWf6RNUNzfOOupswHN9onIDlbSc1N3flfiAXUer95g671d+CJ+bgayQOhtn3U9Ki1WrDdhPVItX3T73tem+dr4Czdelezp1xgMW3wfsOuq8X/khfG4GsgXSiY1zsX5S6pROf6LyTaieVv3kV64sXp38CjS3l/o6sc20c8ZjKerkaU3ozfXTQf0LoFM3aoV6V0YuVZ0cKFr3Cp9ev4igTj/5lSuLV6e+As0r4hrq7A86tc10++KrEjpxtSb07vpZLJBFxOkRMRoR10XEOaXqWAhunPV08v43vgk1dOIO353+yhV1Tqe+As3vra2/P+jUNtPtKyMXo7rve736Ya/IKcuI6APeDzwX+CHwjYj4fGZ+t0Q9nebpyPo6MVAU5v8mtBxO5XTq++3q9vlyPwWzGNXp8zpteuF7a1ttx3X3B53cZpbi2LBOqvu+16sf9kodITsJuC4zr8/M+4FPAWcUqqXjPB3ZOXU/UdXp805/5dNiO3LQqTt8d/IrV9RZnfoKtE5/b223tdr2OvX1Q9DZbUb13vfqrp/LTWRm9xca8VLg9Mz83er3VwG/mJlnz9T+5wcG8t9e+KJuljir23bdx823T3DfnkkOWtHH0av6OfzQg/Zpc/1td7N374N9e8ABweMOP2Sftt+95U4ATjjy4bMus06bTs7ru7fcyZ49ezjx6FVdW95cbW7bdR8/+MndZOa8+nz7TePct+ehby4AB63oY90xA7Xb1F3eVLs660urNnXaXXb9T2fsP4CTH/fItutu1efttHM971ybutvDXG3qrgd11qmp+XVrPV/M23qntoW67Zbret7J/etCOOj4J/Hot71tzjYRcWVmPr2d+ZYKZL8GrJ8WyE7KzN9vanMWcBbA4w95+C989tRnzTnPm+7aC8Axh81+0G++be64P/nR3Xtp7rIIePQhB/CIA2Oftrfdk+zem6w8IDj8YbFPm8VscnKSvr6+1g0XkVZ9Xuf1+97t++6gpzxp1YP98YPxvezeu++2s/KA4PEDB9ReXt11qk67OjXV6adeshTX806psx4s1fW8znbs/nxx68T+fEqd9/667W66ay9jq9bwpD94xZzzOvXUU9sOZKVue/FD4Oim348CdjY3yMzzgfMBDjryuPy9Z75pzrE87/7QpcDs5+W3bh9red66VZtnnnfJjOMEBgf6+fo5p7X6m5eUkZERhoeHS5fRcTONK/nlptf492q+xhvO+QIzfZQJ4Ibzng/UW1/qrlN12t2wfYw/m+E+QOeeuZaneDp8Rst1Pe+UOutUt9fzOtte3e241f5guViO63nddarOe3877aayxuuHOz8GsFQg+wZwXEQcC4wBLwPmjJvzGUw6NZ7g/sm9s86rTpteHWi4nGxYNzjn+lP3ZrydGjjdySsap/6upXhBghanOutUt9fzOtteO18/5PaxNLUzJniu9/V22i20IoP6M3MPcDawDbgWuDAzv9PqeXNdlr39pnEuv+H2/b7Uv1MDorW0TV2aPjjQTzD7pemdGjhdd51q50KRr59zGjec93wvFFFHtFqnur2e19n26m7HWrrqrCudvqdiq6wxX8XuQ5aZ/5aZT8zMx2fmu+o+b7bLsqcn23avpqnTxitpekOdUFNnh9+pK+LaaSd1W7fX87phyw8ny1uddaWTR2brZI35WnJfndTOZdntfGVFnTaeElKzVqc76qwvzW3GxidmHbvguqfFqsR67qlG1VlXOnlPxW58D+eSCmQzfaKqe2Sr1XgCxxxoIdRZX6batBp467qnxcr1XCV0akxwnXbdGENe5LYX7ep72COy79BV90/uun1s78SdtzdPW7n6sWKJuuAAAAWlSURBVGujb8WB05+Tk3vu3/2TG3dM/X5A/8NX9R26ajD6VhyYk3tmnFedNj3kcOC20kX0GPu8++zz7rPPu69n+7zu+3qrdnWzRpOhzDysnVqXxBGyvRN3Xjl5zx1t3c9D8xMR32z3HiqaH/u8++zz7rPPu88+776I+Ga7zyk2qF+SJEkNBjJJkqTClkogO790AT3IPu8++7z77PPus8+7zz7vvrb7fEkM6pckSVrOlsoRMkmSpGVrUQeyiPiLiLg6Iq6KiC9HxJrq8YiI90bEddX0p5WudbmIiM0R8b2qXy+KiIGmaZuqPh+NiPUl61xOIuLXIuI7EbE3Ip4+bZp9vkAi4vSqX6+LiHNK17McRcRHIuLWiLim6bFVEXFxRHy/+vfnSta43ETE0RHxlYi4ttqvvLF63H5fIBFxcERcERHfrvr8ndXjx0bE5VWfXxAR+9w2o9miDmTA5sw8MTOfCvwr8CfV488Djqt+zgI+WKi+5ehi4MmZeSLwX8AmgIg4gcaXwP88cDrwgYjom3Uuasc1wJnAV5sftM8XTtWP76exLzkBeHnV3+qsj9FYd5udA/xHZh4H/Ef1uzpnD/DmzDweOBl4Q7Vu2+8L5z7gtMx8CvBU4PSIOBn4S+D/Vn3+M+C1c81kUQeyzLyz6ddDgKkBb2cA/y8bLgMGIuLIrhe4DGXml6svfwe4DDiq+v8ZwKcy877MvAG4DjipRI3LTWZem5mjM0yyzxfOScB1mXl9Zt4PfIpGf6uDMvOrwPSbcJ4B/GP1/38ENnS1qGUuM2/JzG9V/78LuBYYxH5fMFUW2VX9urL6SeA04DPV4y37fFEHMoCIeFdE3Ay8kgePkA0CNzc1+2H1mDrrd4AvVv+3z7vPPl849m05R2TmLdAID8CjCtezbEXEY4F1wOXY7wsqIvoi4irgVhpnmn4AjDcd4Gi5jykeyCLi3yPimhl+zgDIzLdn5tHAx4Gzp542w6y8XLSmVn1etXk7jUPfH596aIZZ2ec11enzmZ42w2P2eWfYt1rWIuJQ4LPAm6adbdICyMzJanjVUTSOwB8/U7O55lH8q5My8zk1m34C+ALwpzSS5tFN044Cdna4tGWrVZ9HxGuAFwDPzgfvi2Kfz0Mb63kz+3zh2Lfl/DgijszMW6qhJreWLmi5iYiVNMLYxzNzS/Ww/d4FmTkeESM0xu8NRMSK6ihZy31M8SNkc4mI45p+fRHwver/nwdeXV1teTJwx9ShWM1PRJwOvBV4UWbe0zTp88DLIuKgiDiWxgUVV5SosYfY5wvnG8Bx1VVQB9K4eOLzhWvqFZ8HXlP9/zXA5wrWsuxERAAfBq7NzL9pmmS/L5CIWD11R4KI6AeeQ2Ps3leAl1bNWvb5or4xbER8FhgC9gL/Dbw+M8eqFe59NK7euQf47cxs+4s8ta+IuA44CPhp9dBlmfn6atrbaYwr20PjMPgXZ56L2hERLwb+FlgNjANXZeb6app9vkAi4leB9wB9wEcy812FS1p2IuKTwDBwOPBjGmc4tgIXAscANwG/lpnTB/5rP0XELwH/Ceyg8d4J8DYa48js9wUQESfSGLTfR+NA14WZ+ecR8TgaFwytArYDv5mZ9806n8UcyCRJknrBoj5lKUmS1AsMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5mknhMRH42IF0TEQER8sbo5ryQVYyCT1IvW0vhWhM8Bf5GZFxWuR1KP8079knpKRBwA3EXj68Hen5l/WbgkSfIImaSecxywE/gt4PURsbJsOZJkIJPUe9YCF2fmJcA1wKsL1yNJBjJJPWctjSAG8G5gU0SsKFiPJDmGTJIkqTSPkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIK+/87Fl5bxfFXsAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "L = 10000  # number of samples\n",
    "K = 30  # limit for lags in ACF\n",
    "\n",
    "# generate input signal (white Gaussian noise)\n",
    "np.random.seed(2)\n",
    "x = np.random.normal(size=L)\n",
    "# compute system response\n",
    "y = np.convolve(x, [1, 1, 1, 1, 1], mode='full')\n",
    "\n",
    "# compute and truncate ACF\n",
    "acf = 1/len(y) * np.correlate(y, y, mode='full')\n",
    "acf = acf[len(y)-K-1:len(y)+K-1]\n",
    "kappa = np.arange(-K, K)\n",
    "\n",
    "# plot ACF\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.stem(kappa, acf, use_line_collection=True)\n",
    "plt.title('Estimated ACF of output signal $y[k]$')\n",
    "plt.ylabel(r'$\\hat{\\varphi}_{yy}[\\kappa]$')\n",
    "plt.xlabel(r'$\\kappa$')\n",
    "plt.axis([-K, K, 1.2*min(acf), 1.1*max(acf)])\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Derive the theoretic result for $\\varphi_{yy}[\\kappa]$ by calculating the filter-ACF $\\varphi_{hh}[\\kappa]$.\n",
    "* Why is the estimated ACF $\\hat{\\varphi}_{yy}[\\kappa]$ of the output signal not exactly equal to its theoretic result $\\varphi_{yy}[\\kappa]$?\n",
    "* Change the number of samples `L` and rerun the example. What changes?\n",
    "\n",
    "Solution: The filter-ACF is given by $\\varphi_{hh}[\\kappa] = \\text{rect}_N[\\kappa] * \\text{rect}_N[-\\kappa]$. The convolution of two rectangular signals $\\text{rect}_N[\\kappa]$ results in a triangular signal. Taking the time reversal into account yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{hh}[\\kappa] = \\begin{cases} \n",
    "N - |\\kappa| & \\text{for } -N < \\kappa \\leq N \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "\n",
    "for even $N$. The estimated ACF $\\hat{\\varphi}_{yy}[\\kappa]$ differs from its theoretic value due to the statistical uncertainties when using random signals of finite length. Increasing its length `L` lowers the statistical uncertainties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-Correlation Function\n",
    "\n",
    "The cross-correlation functions (CCFs) $\\varphi_{xy}[\\kappa]$ and $\\varphi_{yx}[\\kappa]$ between the in- and output signal of an LTI system $y[k] = \\mathcal{H} \\{ x[k] \\}$ are derived. As for the ACF it is assumed that the input signal originates from a wide-sense stationary real-valued random process and that the LTI system's impulse response is real-valued, i.e. $h[k] \\in \\mathbb{R}$.\n",
    "\n",
    "Introducing the convolution into the definition of the CCF and rearranging the terms yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\varphi_{xy}[\\kappa] &= E \\{ x[k+\\kappa] \\cdot y[k] \\} \\\\\n",
    "&= E \\left\\{ x[k+\\kappa] \\cdot \\sum_{\\mu = -\\infty}^{\\infty} h[\\mu] \\; x[k-\\mu] \\right\\} \\\\\n",
    "&= \\sum_{\\mu = -\\infty}^{\\infty} h[\\mu] \\cdot E \\{ x[k+\\kappa] \\cdot x[k-\\mu] \\} \\\\\n",
    "&= h[-\\kappa] * \\varphi_{xx}[\\kappa]\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "The CCF $\\varphi_{xy}[\\kappa]$ between in- and output is given as the time-reversed impulse response of the system convolved with the ACF of the input signal. \n",
    "\n",
    "The CCF between out- and input is yielded by taking the symmetry relations of the CCF and ACF into account\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{yx}[\\kappa] = \\varphi_{xy}[-\\kappa] = h[\\kappa] * \\varphi_{xx}[\\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "The CCF $\\varphi_{yx}[\\kappa]$ between out- and input is given as the impulse response of the system convolved with the ACF of the input signal. \n",
    "\n",
    "For a system which just attenuates the input signal $y[k] = A \\cdot x[k]$, the CCFs between input and output are given as $\\varphi_{xy}[\\kappa] = A \\cdot \\varphi_{xx}[\\kappa]$ and $\\varphi_{yx}[\\kappa] = A \\cdot \\varphi_{xx}[\\kappa]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System Identification by Cross-Correlation\n",
    "\n",
    "The process of determining the impulse response or transfer function of a system is referred to as *system identification*. The CCFs of an LTI system play an important role in the estimation of the impulse response $h[k]$ of an unknown system. This is illustrated in the following.\n",
    "\n",
    "The basic idea is to use a specific measurement signal as input signal to the system. Let's assume that the unknown LTI system is excited by [white noise](../random_signals/white_noise.ipynb). The ACF of the wide-sense stationary input signal $x[k]$ is then given as $\\varphi_{xx}[\\kappa] = N_0 \\cdot \\delta[\\kappa]$. According to the relation derived above, the CCF between out- and input for this special choice of the input signal becomes\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{yx}[\\kappa] = h[\\kappa] * N_0 \\cdot \\delta[\\kappa] = N_0 \\cdot h[\\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "For white noise as input signal $x[k]$, the impulse response of an LTI system can be estimated by estimating the CCF between its out- and input signals. Using noise as measurement signal instead of a Dirac impulse is beneficial since its [crest factor](https://en.wikipedia.org/wiki/Crest_factor) is limited."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "The application of the CCF to the identification of a system is demonstrated. The system is excited by wide-sense ergodic normal distributed white noise with $N_0 = 1$. The ACF of the in- and output, as well as the CCF between out- and input is estimated and plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\hat{h}[k]$, $h[k]$')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAADgCAYAAABVT4GsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5QkZX3/8fdndwFR1FVZE9gFQUUjEZW4ohxMQoIGJArGGIV4gXhBE40mJipookSjknBO1CR4IWoSlajECxKDwXtM/ImyuiiCogQQdkFugqAgsuz390fXYG9v92z3TM901877dc6c6ap6uupb9dTl20/V052qQpIkSe2wbNIBSJIkaXgmb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIiZvkiRJLWLyJmloSX41yUWTjqOfJAcn2TDpOLYlHf+c5IYkX+0z/RlJPjWJ2IaVZM8kP06yfBGWVUkeuNDLkdrE5E1aApJcluTW5oI78/ePQ7xviwtnVf1PVT14gWL8lyR/vRDz7lpGklyS5MIB0w9N8sUkNye5Nsl/JzmimXZskjtG3YZ9PBZ4PLCmqg7onVhVp1XVb81hviNJsldTvytGfW9VXV5Vu1TVHQsRm6TZjXzQSmqtJ1XVZyYdxIT9GnBfYEWSR1XVuTMTkjwVeA/wMuBJwM3ArwLPBM5sin25qh47zxjuB1xWVT+Z53wkLVG2vElLXJIHNi1MP0pyXZIPNeO/2BT5RtPK9PTeW5NNi97Lk3wzyU+SvDvJLyT5ZNN69Zkk9+oq/+9JftAs64tJfrkZfxzwDOAVzbL+oxm/e5KPNK1glyZ5Sde8dm5a625oWtIeNcTqHgN8HDireT0zrwB/B7y+qt5VVT+qqs1V9d9V9fw5bNPdk5yZ5IdJLk7y/Gb8c4F3AQc26/lXfd57bJL/7RquJC9M8r1mXU9p4p0p+6Uk/9Bs0+8kOaTrvZcleVzX8IlJ3t8MztTvjU0sB/aJ5YAk65LclOTqJH/XjN+i1S7J3l0tlp9pYnx/T9ljklze7GOv7lnGl5PcmOSqJP+YZMdRt7m0lJi8SXo98CngXsAa4B8AqurXmukPb26RfWjA+3+Xzm3AB9Fpsfok8CpgVzrnmJd0lf0ksA+d1q+vA6c1yzq1ef23zbKelGQZ8B/AN4DVwCHAnyQ5tJnXa4EHNH+H0pWM9ZPkrsBTm+WcBhzVlSQ8GNgD+PBs8xjBB4ANwO7NMt+Y5JCqejfwQjoteLtU1WuHnN8T6SSnDweeRmd9ZzwauITO9n4t8NEk9x5injP1u7KJ5ct9yrwVeGtV3YPOdj59wLz+DfgqcB/gROBZfco8ls52PgR4TZKHNOPvAP60if/AZvofDRG/tGSZvElLxxlN68bM30yL0u10buXtXlU/rar/nWUe/fxDVV1dVRuB/wG+UlXrq+o24GPA/jMFq+o9VXVzM+1E4OFJ7jlgvo8CVlXV66rqZ1V1CfBPwFHN9KcBb6iqH1bVFcDfbyPOpwC30UlUP0HnsZHfbqbdp/l/1Tbm8ZiebfiY3gJJ9qCTqLyy2Z7n0Wlt65fQDOukqrqxqi4HPg88omvaNcBbqur2JsG+qGu95ut24IFJdq2qH1fVOb0FkuxJp65e09TT//Lz28zd/qqqbq2qb9BJyB8OUFVfq6pzqmpTVV0GvBP49THFL22XTN6kpePJVbWy6++fmvGvAAJ8NckFSZ4z4nyv7np9a5/hXQCSLE9yUpL/S3ITcFlTZtcB870fsHt3skSnRe8Xmum7A1d0lf/+NuI8Bji9SRJuAz7Kz1vrrm/+77aNeZzTsw23SmaauH5YVTf3xLZ6G/OezQ+6Xt9Cs00bG6uqepa1+zyW1e25dFpUv5Pk3CRP7FNmZn1v6Rp3RZ9yfdchyYOSfKK5nX4T8EYG7xOSMHmTlryq+kFVPb+qdgdeALwtC/PVDL8PHAk8DrgnsFczPjOh9JS/Ari0J1m6e1Ud3ky/is6tzhl7DlpwkjXAbwLPbJKEH9C5nXl4kl3ptFZdQecW8HxdCdw7yd17Yts4hnn3s3rmGbiuZV3ZvP4JcNeuab/Y9bp3e2+lqr5XVUfTuc39N8CHk9ytp9hVdNa3ezl7MLy3A98B9mluz76Kn+8TkvoweZOWuCS/1yQ3ADfQuajPfAXE1cD9x7Sou9O5bXk9nYTijT3Te5f1VeCmJK9sOicsT/LQJDMdE04HTkhyryb+P55l2c8CvkvnmatHNH8PovNc2tFNy9XLgL9M8gdJ7pFkWZLHJjl1lJVsbuH+P+BNSe6S5GF0WrBOG2U+I7gv8JIkOyT5PeAhdDpkAJxH59m+HZKspZOwzrgW2Mws9ZvkmUlWVdVm4MZm9BZfD1JV3wfWAScm2bHp+PCkEeK/O3AT8OMkvwT84QjvlZYkkzdp6fiPbPkdZR9rxj8K+EqSH9N5VumlVXVpM+1E4F+b25ZPm+fy30vnlt5G4EKg95bju4F9m2Wd0XyH2JPoJFqXAtfReXZs5hm5v2rmdymd59jeN8uyjwHe1rQy3vkHvKOZRlV9GHg68Bw6LVdXA39Np3fqqI6m07J4JZ3n/l5bVZ+ew3yG8RU6nUCuA94APLWqZm4D/yWdjgY30Nle/zbzpuY25xuALw16fg84DLig2TfeChxVVT/tU+4ZdDobXE9nm32ITqI+jD+n0yp7M51nGgd1jJHUyJaPSkiS2iLJscDzxvDdc2OVztfNfGeE3rSSRmDLmyRpXpI8KskDmlvNh9F5tvGMScclba/8hQVJ0nz9Ip3eu/eh8xzhH1bV+smGJG2/vG0qSZLUIt42lSRJahGTN0mSpBZZUs+87brrrrXXXntNOgxJkqRt+trXvnZdVa3qHb+kkre99tqLdevWTToMSZKkbUrS92f/vG0qSZLUIiZvkiRJLWLyJkmS1CImb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIlOZvCV5T5JrknxrwPQk+fskFyf5ZpJfWewYJemM9Rs56KTPsffx/8lBJ32OM9ZvnHRIkpaAqUzegH8BDptl+hOAfZq/44C3L0JMknSnM9Zv5ISPns/GG2+lgI033soJHz3fBE7SgpvK5K2qvgj8cJYiRwLvrY5zgJVJdluc6CQJTj77Im69/Y4txt16+x2cfPZFE4pI0lIxlcnbEFYDV3QNb2jGbSXJcUnWJVl37bXXLkpwkrZ/V95460jjJWlc2pq8pc+46lewqk6tqrVVtXbVqq1+21WS5mT3lTuPNF6SxqWtydsGYI+u4TXAlROKRdIS9PJDH8zOOyzfYtzOOyzn5Yc+eEIRSVoq2pq8nQk8u+l1+hjgR1V11aSDkrR0PHn/1bzpKfux4/LOaXT1yp1501P248n7932CQ5LGZsWkA+gnyQeAg4Fdk2wAXgvsAFBV7wDOAg4HLgZuAf5gMpFKWsqevP9qPvDVywH40AsOnHA0kpaKqUzequrobUwv4EWLFI4kSdLUaOttU0mSpCXJ5E2SJKlFTN4kSZJaxORNkiSpRUzeJEmSWsTkTZIkqUVM3iRJklrE5E2SJKlFTN4kSZJaxORNkiSpRUzeJEmSWsTkTZIkqUVM3iRJklrE5E2SJKlFTN4kSZJaxORNkiSpRaY2eUtyWJKLklyc5Pg+0/dM8vkk65N8M8nhk4hTkiRpMU1l8pZkOXAK8ARgX+DoJPv2FPsL4PSq2h84Cnjb4kYpSZK0+KYyeQMOAC6uqkuq6mfAB4Eje8oUcI/m9T2BKxcxPkmSpIlYMekABlgNXNE1vAF4dE+ZE4FPJflj4G7A4xYnNEmSpMmZ1pa39BlXPcNHA/9SVWuAw4H3JdlqfZIcl2RdknXXXnvtAoQqSZK0eKY1edsA7NE1vIatb4s+FzgdoKq+DNwF2LV3RlV1alWtraq1q1atWqBwJUmSFse0Jm/nAvsk2TvJjnQ6JJzZU+Zy4BCAJA+hk7zZtCZJkrZrU5m8VdUm4MXA2cC36fQqvSDJ65Ic0RT7M+D5Sb4BfAA4tqp6b61KkiRtV6a1wwJVdRZwVs+413S9vhA4aLHjkiRJmqSpbHmTJElSfyZvkiRJLWLyJkmS1CImb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIiZvkiRJLWLyJkmS1CImb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIiZvkiRJLWLyJkmS1CImb5IkSS0ylclbksOSXJTk4iTHDyjztCQXJrkgyb8tdoySJEmTsGLSAfRKshw4BXg8sAE4N8mZVXVhV5l9gBOAg6rqhiT3nUy0kiRJi2saW94OAC6uqkuq6mfAB4Eje8o8Hzilqm4AqKprFjlGSZKkiZjG5G01cEXX8IZmXLcHAQ9K8qUk5yQ5bNGikyRJmqCpu20KpM+46hleAewDHAysAf4nyUOr6satZpYcBxwHsOeee443UkmSpEU2jS1vG4A9uobXAFf2KfPxqrq9qi4FLqKTzG2lqk6tqrVVtXbVqlULErAkSdJimcbk7VxgnyR7J9kROAo4s6fMGcBvACTZlc5t1EsWNUpJkqQJmLrkrao2AS8Gzga+DZxeVRckeV2SI5piZwPXJ7kQ+Dzw8qq6fjIRS5IkLZ5pfOaNqjoLOKtn3Gu6XhfwsuZPkiRpyRgqeUty7yGKbe7XYUCSJEnjM2zL25XNX7+eoDOWA3bnlCRJWkDDJm/frqr9ZyuQZP0Y4pEkSdIshu2wcOCYykiSJGkehkrequqnAEl6f6aKJMu6y0iSJGnhjPpVIS9I8mjo/IB8kucA3xl/WJIkSepn1K8KORr4eJL/BP4QOB949tijkiRJUl+jJm8HAa8GTgOOraovjD0iSZIkDTRq8vZUYD9gV+CdST4NfKuq3jH2yCRJkrSVkZK3qnoOQJLQ+SH4/Zo/SZIkLYI5/TxW8/NU323+PjLWiCRJkjTQUL1Nk3x9HGUkSZI0P8O2vD0kyTdnmR7gnmOIR5IkSbMYNnn7pSHK3DGfQCRJkrRtQyVvVfX9fuOTLK8qkzZJkqRFMuovLPQ6NcldAZL82hjikSRJ0izmm7y9Bnh3kvcBjxpDPHdKcliSi5JcnOT4Wco9NUklWTvO5UuSJE2j+SZvrwcuAgo4ff7hdCRZDpwCPAHYFzg6yb59yt0deAnwlXEtW5IkaZrNN3l7RVWdSOd3Tl87/3DudABwcVVdUlU/Az4IHNmn3OuBvwV+OsZlS5IkTa35Jm8nJdm5qn4CfGAcATVWA1d0DW9oxt0pyf7AHlX1iTEuV5IkaarN6RcWurwWeE+STcB5wGfnHxLQ+d64XnXnxGQZ8Gbg2G3OKDkOOA5gzz33HFN4kiRJkzGVz7zRaWnbo2t4DXBl1/DdgYcCX0hyGfAY4Mx+nRaq6tSqWltVa1etWjXGECVJkhbfyC1vSR4F7AhcQOeZt+uS3A14K/C8McV1LrBPkr2BjcBRwO/PTKyqHwG7dsX0BeDPq2rdmJYvSZI0lUZK3pKcQacV7Ao6vUC/lOTFVfWTJC8YV1BVtSnJi4GzgeXAe6rqgiSvA9ZV1ZnjWpYkSVKbzJq8NV/P8aqqemYz6qHAM4D1VfWzJM8B3g4cM+5fWqiqs4Czesa9ZkDZg8e5bEmSpGm1rZa3zwIHdg2fBLwSeHiSW4DzgV9P8njgvKq6dmHClCRJEmy7w8JvAW/oGn4ScHpVPQB4LPAOOrc1jwL+a0EilCRJ0p1mbXmrqvPp3Cad8Tzgn5OcAHwLeDDwX1X13IULUZIkSTNG6rDQ3BZ9YpLdgf2AH1XVOQsSmSRJkrYypy/praor2fJ71yRJkrQI5vslvZIkSVpEJm+SJEktYvImSZLUIiZvkiRJLWLyJkmS1CImb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIiZvkiRJLWLyJkmS1CJTmbwlOSzJRUkuTnJ8n+kvS3Jhkm8m+WyS+00iTkmSpMU2dclbkuXAKcATgH2Bo5Ps21NsPbC2qh4GfBj428WNUpIkaTKmLnkDDgAurqpLqupnwAeBI7sLVNXnq+qWZvAcYM0ixyhJkjQR05i8rQau6Bre0Iwb5LnAJxc0IkmSpCmxYtIB9JE+46pvweSZwFrg1wfOLDkOOA5gzz33HEd8kiRJEzONLW8bgD26htcAV/YWSvI44NXAEVV126CZVdWpVbW2qtauWrVq7MFKkiQtpmlM3s4F9kmyd5IdgaOAM7sLJNkfeCedxO2aCcQoSZI0EVOXvFXVJuDFwNnAt4HTq+qCJK9LckRT7GRgF+Dfk5yX5MwBs5MkSdquTOMzb1TVWcBZPeNe0/X6cYselCRJ0hSYupY3SZIkDWbyJkmS1CImb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIiZvkiRJLWLyJkmS1CImb5IkSS1i8iZJktQiJm+SJEktYvImSZLUIiZvkiRJLbJi0gFIS9kZ6zdy8tkXceWNt7L7yp15+aEP5sn7r550WNJ2z2NPbWbyJk3IGes3csJHz+fW2+8AYOONt3LCR8+/c7oXFmk8ehO13/ilVXzkaxv7HnseZ2oDk7cefhr7uVG2xVLcbvNd55PPvujOi8eMW2+/gxPPvIDbNm3e6sKy7vs/5PPfuXar5fWLY2b+w5SdbfxCrPcgS3EfWmwLuY3nuw/B/PfZfvMAtvqQdNo5l1M9Md16+x2cfPZFI22PpbjPjrrOS3EbDTLObZGq3l14OiQ5DHgrsBx4V1Wd1DN9J+C9wCOB64GnV9Vls81zp932qbUvfefQBznAzjss501P2a/vBp6GC94g811eb6sQDN4Ws5WFhWlBGsf2HOUC0u+9o+wr/ex9/H9udQGZTWCL8jvvsJzffeTqLVoQAHZYFgjcfkdts+xs40et62HXe9B2X6h9aKH2lZl5PP2dXwbgQy84cFGWN9d59LY2wXi38XzOF+PYZwfN4y47LOOGW24faj0CXHrSb499nWfKz/eD02yxLMZ5drZ9aKHOk9N87Ry1bL9t8buPXN33Q/mMJF+rqrW985vK5C3JcuC7wOOBDcC5wNFVdWFXmT8CHlZVL0xyFPA7VfX02ea702771G7HvGXkg3z1yp350vG/ucW4cSU3i32AQf+ktXfcyWdfxMYbb91qWSt33oG77bRi6LLdLUjbWu+5fpKeme+2DoLeZQ17AekX70Enfa7vOvfbVwYZNI9RLE+4Y8hjeFDZQeNHqeth13vQsTDo2Bu0Dw2q63EkLKNesLqTt1E/yAwatxDz6E3+t7WNRzk/DdqXR9mH+hl1n52vcRy/o1wvRvngNDOf+Z4Ph004+sU8aB9a3VW3c72OjOv4nY+FvK4P2l/6fSjvnkfbkrcDgROr6tBm+ASAqnpTV5mzmzJfTrIC+AGwqmZZoZnkbS5Wr9x56ItY7058y882DUwK++3wo95qmM9JelDC0ns7b5BRys4Y9uQ2apI96CCA4ZPTQfH21tOg9872yX3YE9MoLQWLaVt13XuMjHICG9WwrZCjJiyjzGNmP+5O3mZLYoY99kZJZMfR2jTIoESm33noTz903lCtyHM5XyyUUc4X/fblQS3nAd789EcMdQ0Y5YMT9E/SRjkfjpIsjnqs9tbtqNeRcR2/w36Ah6335VGu1aNe14c9RmbmMXPstS15eypwWFU9rxl+FvDoqnpxV5lvNWU2NMP/15S5btB8H7DLPeuND3/svONbtixs3jx4u21r+mxlly0Lq3bZiRtvuZ3bNt3BTiuWs8e9dwbgkut+slXZZQmb7tg8xzUZLAnD7hujlJ2x04rlW6zfFT+8lds2jf+kvmL5MjZXbbXdhq2fUd/Tuz7bqr9Vu+zENTffRlXNWnaQcdTTQtT1smXh/rveDWCLbbEQdTxqbOOcx2Pufx8uvOomAPbd7R6cc8n181r+tHngfXfZov5W3nUHrv3xbfM6Dy32PjvoHDDssTfonDzonNVvefM1rnP9oO22YvkylidbrN/F1/x43vNd6OvItgw6D416np2Z13yu69uqv0vuuZp3PuxIYMtGgEHJ27R2WEifcb1bbZgyJDkOOA7g/ne7x6wLHfag27y5Zt1Zh63gfmU3by6uvumndw7ftukOLrnuJywbUHbzSE9NDa+qht5Z+5Xd1s46c9KbWb9xnui69Vv+bPXXz7B1umxZWHnXHbZYn23V34233M6j97533/kNc9FctctOW41POodG9/oNKjto/CDD7hebNxeXXX/LFsfTbInboAvsKBescZz4R5nHTiuWA52krXvcQiao4zbbNl6xfNlW+/LVN229bps3F8uWZ17ni3Hss4Pmsdd97gpsffHedZed2Lu5sM9Yf/mNQ5+TOzFs3iq2mfcMY9jz0LjO9YOWtemOzWxqXs+s34rly4Y69mar91GvI+M26Dw06Jw8yDiu66McI7uv3HmbZaY1edsA7NE1vAa4ckCZDc1t03sCP+ydUVWdCpwKndumr/zVP5r12ablMNTtsZn3zKeZeCFvH/RrMh/lVso4mokB/mLIW7qjPMPSr8l80Hxn01sHc7mF3Hur8NUj3I6dibvfbdb70emJ063f7aqDR+h1N6hsv/Gj1PWot0EH3a7qPfbGtQ+NciyMMo83PWU/7tdzS2b9+o2cOMKtrX4G3RIaxzz63VaCrbfxXB70H/ZW4Wzni/nss4Pm8cimjnqPp34OHqET0Xxvjw26VTiqUc6H8z3PDtqHxvEY0TiO34Uwrut6v2Nk0KMzM/vybKb1tukKOh0WDgE20umw8PtVdUFXmRcB+3V1WHhKVT1ttvluq7fpKM/ojOMBzVGeuxpk1JP0fB7Sn0vHi94T7CjJ8GyxAfN6fmyUC8goD+mP2oN0lAekF9ModT2X59iGeT6uO5Zh6nq2nonDHAujzmO+z6aOsn+Pax6jxDzXZ3S657mQHbUWwij78qAPXqN03uj3kP6g5Ga2jiUw3PlwlGSxX7Ixjt73o3beGPb4ncsH+H5G7ZQ3305d2+pA0qpn3gCSHA68hc5Xhbynqt6Q5HXAuqo6M8ldgPcB+9NpcTuqqi6ZbZ5r166tdevWjRTHXLqDz2cnnkuPMJjfSXq+7x/lRDxqMjyOdRtHV/X5JjLj6NG32ObTK23U3tvjim0c3fwX83vsYLQec+OYx7Dm2jtuW/FO6/4Oo/ewHHQxns85Z1xfvTTs/j1bC+kox+liH3vz7QA2yjl5kt8u0brkbSHMJXmDhf1Swkl3jV5sk/g0vpjfozSJ77ybBguVOGty5vq9VG033+826zePSXw34SjL2l6O04X8DslJfa+ryRtzT94WW9s+rY7K9Vs63BbtZv11bO/bwfWbXiZvtCd5kyRJGpS8LZtEMJIkSZobkzdJkqQWMXmTJElqEZM3SZKkFjF5kyRJahGTN0mSpBYxeZMkSWoRkzdJkqQWMXmTJElqEZM3SZKkFjF5kyRJahGTN0mSpBYxeZMkSWoRkzdJkqQWmbrkLcm9k3w6yfea//fqU+YRSb6c5IIk30zy9EnEKkmStNimLnkDjgc+W1X7AJ9thnvdAjy7qn4ZOAx4S5KVixijJEnSRExj8nYk8K/N638FntxboKq+W1Xfa15fCVwDrFq0CCVJkiZkGpO3X6iqqwCa//edrXCSA4Adgf9bhNgkSZImasUkFprkM8Av9pn06hHnsxvwPuCYqto8oMxxwHHN4I+TXDTKMlpmV+C6SQehObHu2s36ay/rrt229/q7X7+RqarFDmRWTXJ1cFVd1SRnX6iqB/cpdw/gC8CbqurfFznMqZRkXVWtnXQcGp11127WX3tZd+22VOtvGm+bngkc07w+Bvh4b4EkOwIfA95r4iZJkpaSaUzeTgIen+R7wOObYZKsTfKupszTgF8Djk1yXvP3iMmEK0mStHgm8szbbKrqeuCQPuPXAc9rXr8feP8ih9YGp046AM2Zdddu1l97WXfttiTrb+qeeZMkSdJg03jbVJIkSQOYvG0nkvx5kkqyazOcJH+f5OLmJ8R+ZdIxamtJTk7ynaaOPtb9SyFJTmjq76Ikh04yTvWX5LCmfi5O0u/XYDRFkuyR5PNJvt38vOJLm/Hb/FlGTYcky5OsT/KJZnjvJF9p6u5DTYfG7Z7J23YgyR50Ondc3jX6CcA+zd9xwNsnEJq27dPAQ6vqYcB3gRMAkuwLHAXM/ATc25Isn1iU2kpTH6fQOdb2BY5u6k3TaxPwZ1X1EOAxwIuaOhvmZxk1HV4KfLtr+G+ANzd1dwPw3IlEtchM3rYPbwZeAXQ/wHgkna9Sqao6B1jZfG+epkhVfaqqNjWD5wBrmtdHAh+sqtuq6lLgYuCAScSogQ4ALq6qS6rqZ8AH6dSbplRVXVVVX29e30wnCVjNED/LqMlLsgb4beBdzXCA3wQ+3BRZMnVn8tZySY4ANlbVN3omrQau6Bre0IzT9HoO8MnmtfU3/ayjFkuyF7A/8BVG/FlGTcxb6DRUzPyi0n2AG7s+AC+ZY3DqvipEW9vGz4m9Cvitfm/rM86uxRMwW/1V1cebMq+mc0vntJm39Slv/U0X66ilkuwCfAT4k6q6qdOAo2mW5InANVX1tSQHz4zuU3RJHIMmby1QVY/rNz7JfsDewDeak88a4OtJDqDzCWSPruJrgCsXOFT1Maj+ZiQ5BngicEj9/Lt7rL/pZx21UJId6CRup1XVR5vRVyfZretnGa+ZXIQa4CDgiCSHA3cB7kGnJW5lkhVN69uSOQa9bdpiVXV+Vd23qvaqqr3oXEx+pap+QOdnxp7d9Dp9DPCjmdsCmh5JDgNeCRxRVbd0TToTOCrJTkn2ptPx5KuTiFEDnQvs0/R225FOB5MzJxyTZtE8I/Vu4NtV9Xddk7b5s4yarKo6oarWNNe6o4DPVdUzgM8DT22KLZm6s+Vt+3UWcDidB91vAf5gsuFogH8EdgI+3bSenlNVL6yqC5KcDlxI53bqi6rqjgnGqR5VtSnJi4GzgeXAe6rqggmHpdkdBDwLOD/Jec24V9H5GcbTkzyXTq/935tQfBrdK4EPJvlrYD2d5Hy75y8sSJIktYi3TSVJklrE5E2SJKlFTN4kSZJaxORNkiSpRUzeJEmSWsTkTZIkqUVM3iRJklrE5E2S5ijJP6TIIFMAAAC8SURBVCd5YpKVST6Z5HcmHZOk7Z/JmyTN3X7AjXR+kuf1VfWxCccjaQnwFxYkaQ6SLANuBq4HTqmqv5lwSJKWCFveJGlu9gGuBI4FXphkh8mGI2mpMHmTpLnZD/h0VX0O+Bbw7AnHI2mJMHmTpLnZj07SBvBG4IQkKyYYj6QlwmfeJEmSWsSWN0mSpBYxeZMkSWoRkzdJkqQWMXmTJElqEZM3SZKkFjF5kyRJahGTN0mSpBYxeZMkSWqR/w9/wCgf2AUeuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAADgCAYAAAC3k5rVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5xcdX3v8dc7myWsRlkwW002IaGabkWjRiKFi+21/mgoCqRIJdYqKDbQwlXvbYMGW4u0Bby51/b6uwhUVC4EMcZooSkWrdIrPwIBQsDYiALZ8CP8WCCyhmTzuX+cM2EyOTM7szs758zO+/l4zGN3vue7Zz7nnJk9n/me7/d7FBGYmZmZWXFMyTsAMzMzM9uXEzQzMzOzgnGCZmZmZlYwTtDMzMzMCsYJmpmZmVnBOEEzMzMzKxgnaGYdStJvS9qcdxxZJL1J0ta84xiNEv8k6UlJt+Qdz0SRdK6kS1rwOqdJunGiX8esHThBM2szkn4haVjSjrLH5+r4u5D0itLziPhRRAxMUIxfkfS3E7HusteQpPsk3VNl+WJJP5T0jKTtkv5d0gnpstMkjTS6DzO8EXgbMDsijhzH5mTFf56krxdhfRFxQUR8sFmxmNnopuYdgJmNyfER8b28g8jZ7wC/BkyV9IaIuLW0QNLJwGXA/wCOB54Bfhv4Y2BtWu3HEfHGccYwF/hFRPxynOsxM9uHW9DMJhFJr0hbip6S9JikVWn5D9Mqd6atRadUXkZMW+aWS7pL0i8lXSrppZKuS1uhvifp4LL635D0cPpaP5T0qrR8GfAe4Jz0tb6Tls+S9M20Nevnkj5Utq6etNXtybRF7A11bO6pwLeBa9PfS+sS8GngbyLikoh4KiL2RMS/R8SfjGGfzpK0VtITkrZI+pO0/HTgEuDodDs/mfG3UyT9paT7JT0q6auSDkqX7XcZNz0Gb5V0LHAucEq67jvT5T+QdKGkW9L9/m1Jh4x1fRnxflTSYHq8N0t6S1q+T+ubpPel2/S4pL8qvU5Z3avTbX1G0iZJi8r+9mOSfpYuu0fSHzR6TMw6gRM0s8nlb4B/BQ4GZgOfBYiI30mXvzYipkfEqip//06SS3a/QdLydB3JiX0Gyf+LD5XVvQ6YT9KKdTtwRfpaF6e//8/0tY6XNAX4DnAn0A+8BfiIpMXpuv4aeHn6WExZwpVF0guAk9PXuQJYKumAdPEAMAe4ptY6GnAlsBWYlb7mBZLeEhGXAmeStMRNj4i/zvjb09LH7wK/DkwHRr2UGhH/AlwArErX/dqyxe8DPpDGsxv4zDjXB4CkAeBs4A0R8SKS4/CLjHqHA18gScJnAgeRHNNyJwBXAb0kLZbl2/wzktbMg4BPAl+XNHO0bTDrNE7QzNrTGklDZY9Sy9AukstusyLiVxHRaIfrz0bEIxExCPwIuDkiNkTETuBbwMJSxYi4LCKeSZedB7y21DqU4Q1AX0ScHxHPRcR9wJeBpenydwF/FxFPRMSDjJ50nATsJElGv0vSXePt6bKXpD8fGmUdR1Xsw6MqK0iaQ9LP7KPp/ryDpNXsvaOsu+Q9wKcj4r6I2AGsIEkmx9O95GsRcXd6WfWvgHdJ6hrH+kpGgGnA4ZK6I+IXEfGzjHonA9+JiBsj4jngE0DlTZ1vjIhrI2IE+BqwNyGMiG9ExLa0VXMV8J9AU/vvmU0GTtDM2tOSiOgte3w5LT8HEHBLemnpAw2u95Gy34cznk8HkNQl6aL0UtXTPN/SMqPKeucCs8oTIpKWuZemy2cBD5bVv3+UOE8Fro6I3WmCuJrnW90eT3+O1ipzU8U+vCmjzizgiYh4piK2yhajamax77bcT5JMvjS7el0q91M31fd73SJiC/ARkmT7UUlXSZqVUXWfYxURz/L8Pi95uOz3Z4EDS0lpenn0jrL3waubEb/ZZOMEzWwSiYiHI+JPImIWcAbwBZWN3GyiPwJOBN5KcqlqXlquUigV9R8Efl6REL0oIo5Llz9Eclmy5NBqLyxpNvBm4I/TPnAPk7TqHCdpBrA5fb13jnnrnrcNOETSiypiG2zg7+dW/O1uksT3l8ALSgvSVrC+srqV+7Ckcj/tAh4bx/qerxDxf9OBE3PT+p/KqPYQyeXz0uv08HyrZU2S5pK0nJ4NvCQieoG7ef59Y2YpJ2hmk4ikP0wTGIAnSU6yI+nzR0j6QTXDi0guMT5OkhRcULG88rVuAZ5OO6H3pC1wr5ZUGgxwNbBC0sFp/P+txmu/F/gpSV+z16WP3yDpJ/buiAiS0Zt/Jen9kl6cdtZ/o6SLG9nI9HLr/wMulHSgpNcAp5P2t6vDlcB/l3SYpOk83w9sd7oNB0p6u6Ru4C9JLjGWPALMS/vvlftjSYen/fDOB65JLyWOdX1A0gdN0pslTQN+RdJiOpJR9RrgeEn/Je3390nqT7BeSPKe3J6+5vtJWtDMrIITNLP29B3tO4fXt9LyNwA3S9pB0jn7wxHx83TZecDl6aWld43z9b9KcnltELgHqLw8eClJX6YhSWvSBOJ4kmTq5yQtPpeQtL5BcpK/P132ryT9lqo5FfhC2lq49wF8KV1GRFwDnELSmX4bSXLytySjPhv1bpIWwm0k/fD+OiKur/NvL0u35Yck2/Yr0uQzIp4C/oxkPwyStICVj8L8RvrzcUm3l5V/DfgKyWXEA0kHboxjfSXTgItIjs3DJIM/zq2sFBGb0m24iqQ17RngUZKEvaaIuAf438CPSY7JAuA/Rvs7s06k5MummZkVnaQfAF+PiAmf1b9eacvgEDC/7MuAmY2TW9DMzKwhko6X9AJJLwT+F7CRjCk5zGzsnKCZmVmjTiS55LuNZC68peHLMWZN5UucZmZmZgXjFjQzMzOzgnGCZmZmZlYw47ndSOHMmDEj5s2bl3cYZmZmZqO67bbbHouIvqxlkypBmzdvHuvXr887DDMzM7NRSap6Wztf4jQzMzMrmFwStPSWKbdIujO9ofMnM+qcJml7elPdOyR9MI9YzczMzFotr0ucO4E3R8SO9J5xN0q6LiIqbxezKiLOziE+MzMzs9zkkqClExruSJ92pw9PyGZmZmZGjn3QJHVJuoPkJrvXR8TNGdXeKekuSddImlNlPcskrZe0fvv27RMas5mZmVkr5JagRcRIRLwOmA0cKenVFVW+A8yLiNcA3wMur7KeiyNiUUQs6uvLHKlqZmZm1lZyH8UZEUPAD4BjK8ofj4id6dMvA0e0ODQzMzOzXOQ1irNPUm/6ew/wVuAnFXVmlj09Abi3dRGamZmZ5SevUZwzgcsldZEkiVdHxHclnQ+sj4i1wIcknQDsBp4ATsspVjMzM7OWUjKgcnJYtGhR+E4CZmY2mazZMMjKdZvZNjTMrN4eli8eYMnC/rzDsiaQdFtELMpaNqlu9WRmZjaZrNkwyIrVGxneNQLA4NAwK1ZvBHCSNsnlPkjAzMzMsq1ct3lvclYyvGuEles25xSRtYoTNDMzs4LaNjTcULlNHk7QzMzMCmpWb09D5TZ5OEEzMzMrqOWLB+jp7tqnrKe7i+WLB3KKyFrFgwTMzMwKqjQQ4Jxr7uK5kT30exRnx3CCZmZmVgDVptNYsrCfK295AIBVZxxds65NHk7QzMzMctbIdBqeeqMzuA+amZlZzhqZTsNTb3QGJ2hmZmY5a2Q6DU+90RmcoJmZmeWskek0PPVGZ3CCZmZmlrNGptPw1BudIZcETdKBkm6RdKekTZI+mVFnmqRVkrZIulnSvNZHamZmNvGWLOznwpMWcEBXclru7+3hwpMWZHb6b6Suta+8RnHuBN4cETskdQM3SrouIm4qq3M68GREvELSUuBTwCl5BGtmZjbRsqbTaEZda0+5tKBFYkf6tDt9REW1E4HL09+vAd4iSS0K0czMzCw3ufVBk9Ql6Q7gUeD6iLi5oko/8CBAROwGngJekrGeZZLWS1q/ffv2iQ7bzMzMbMLllqBFxEhEvA6YDRwp6dUVVbJayypb2YiIiyNiUUQs6uvrm4hQzczMzFoq91GcETEE/AA4tmLRVmAOgKSpwEHAEy0NzszMzCwHeY3i7JPUm/7eA7wV+ElFtbXAqenvJwM3RMR+LWhmZmZmk01eozhnApdL6iJJEq+OiO9KOh9YHxFrgUuBr0naQtJytjSnWM3MzMxaKpcELSLuAhZmlH+i7PdfAX/YyrjMzMzMiiD3PmhmZmZmti8naGZmZmYF4wTNzMzMrGCcoJmZmZkVTF6jOM3MzDrWmg2DrFy3mW1Dw8zq7WH54oGm3Ox8otZrrecEzczMrIXWbBhkxeqNDO8aAWBwaJgVqzdO6HqdpLUfX+I0MzNroZXrNu9NokqGd42wct3mQq7X8uEEzczMrIW2DQ03VJ73ei0fTtDMzMxaaFZvT0Plea/X8uEEzczMrIWWLx6gp7trn7Ke7i6WLx4o5HotH07QzMzMWmjJwn4uPGkBB3Qlp+D+3h4uPGnBuDvyT9R6LR8exWlmZtZiSxb2c+UtDwCw6oyjC79ea71cWtAkzZH0fUn3Stok6cMZdd4k6SlJd6SPT2Sty8zMzGyyyasFbTfw5xFxu6QXAbdJuj4i7qmo96OIeEcO8ZmZmZnlJpcWtIh4KCJuT39/BrgX8EVyMzMzMwowSEDSPGAhcHPG4qMl3SnpOkmvamlgZmZmZjnJdZCApOnAN4GPRMTTFYtvB+ZGxA5JxwFrgPkZ61gGLAM49NBDJzhiMzMzs4mXWwuapG6S5OyKiFhduTwino6IHenv1wLdkmZk1Ls4IhZFxKK+vr4Jj9vMzMxsouU1ilPApcC9EfHpKnVeltZD0pEksT7euijNzMzM8pHXJc5jgPcCGyXdkZadCxwKEBFfAk4G/lTSbmAYWBoRkUewZmZmZq2US4IWETcCGqXO54DPtSYiMzMzs+LIfRSnmZmZme3LCZqZmZlZwThBMzMzMysYJ2hmZmZmBeMEzczMzKxgnKCZmZmZFYwTNDMzM7OCcYJmZmZmVjBO0MzMzMwKJq9bPZmZmU16azYMsnLdZrYNDTOrt4fliwdYsrC/42KwxjlBMzMzmwBrNgyyYvVGhneNADA4NMyK1RsBWpYgFSEGGxtf4jQzM5sAK9dt3psYlQzvGmHlus0dFYONTS4JmqQ5kr4v6V5JmyR9OKOOJH1G0hZJd0l6fR6xmpmZjcW2oeGGyidrDDY2ebWg7Qb+PCJeCRwFnCXp8Io6vw/MTx/LgC+2NkQzM7Oxm9Xb01D5ZI3BxqbuBE3SIXU8eutZV0Q8FBG3p78/A9wLVF4MPxH4aiRuAnolzaw3XjMzszwtXzxAT3fXPmU93V0sXzzQUTHY2DQySGBb+lCNOl3AoY0EIGkesBC4uWJRP/Bg2fOtadlDjazfzMwsD6VO+OdccxfPjeyhP4cRlEWIwcamkQTt3ohYWKuCpA2NvLik6cA3gY9ExNOVizP+JDLWsYzkEiiHHtpQbmhmZjahlizs58pbHgBg1RlHd2wM1rhG+qDVc1TrPvKSukmSsysiYnVGla3AnLLns0la8PYRERdHxKKIWNTX11fvy5uZmZkVVt0JWkT8CkDSiZXLJE0przMaSQIuJWmV+3SVamuB96WjOY8CnooIX940MzOzSW8sE9WeIenhiLhZUhdwKvAx4DcaWMcxwHuBjZLuSMvOJe2/FhFfAq4FjgO2AM8C7x9DrGZmZmZtZywJ2ruBb0v6Z+BPgY3A+xpZQUTcSO3BBkREAGeNIT4zMzOztjaWBO0Y4OPAFcBpEfGDpkZkZmZm1uHGkqCdDCwAZgD/KOl64O70sqSZmZmZjVPDCVpEfAD2dvSfT5KsLWhyXGZmZmYdaywtaMDePmI/TR/fbFpEZmZmZh2ukVs93d6MOmZmZmZWWyMtaK+UdFeN5QIOGmc8ZmZmZh2vkQTtN+uoMzLWQMzMzMwsUXeCFhH3Z5VLmhYRO5sXkpmZmVlna+RenNV8QdLbmrAeMzMzM6MJCVpEnA4cKukzkmY0ISYzMzOzjjbuBE3SYuAw4OXAJZL+YNxRmZmZmXWwZlzinAdcGhFvj4glwJubsE4zMzOzjtWMBO1I4BVlzz8+2h9IukzSo5LurrL8TZKeknRH+vhEE+I0MzMzawtN74MWEU/X8WdfAY4dpc6PIuJ16eP88cZpZmZm1i5y6YMWET8Enhjva5uZmZlNRmO6F6ekNwAHAJuAlwGXRcR96bLPAt9qQmxHS7oT2Ab8RURsasI6zczMmm7NhkFWrtvMtqFhZvX2sHzxAEsW9ucdVk3tGHMnaThBk7QGmA08CBwO/AfwjbIqo/ZBq8PtwNyI2CHpOGANML9KPMuAZQCHHnpoE17azMysfms2DLJi9UaGdyU30xkcGmbF6o0AhU142jHmTjPqJU5Jh0v6elnRq4GzgFMiYgC4EfhiaWGdfdBqioinI2JH+vu1QHe1OdYi4uKIWBQRi/r6+sb70mZmZg1ZuW7z3kSnZHjXCCvXbc4potG1Y8ydpp4WtH8Dji57fhHwUeC1kp4FNgL/Nb2bwB0RsX28QUl6GfBIRISkI0kSycfHu14zM7Nm2zY03FB5EbRjzJ2mnkECvwf8Xdnz44GrI+LlwBuBLwFdwFLgX+p5UUlXAj8GBiRtlXS6pDMlnZlWORm4O+2D9hlgaUREXVtkZmbWQrN6exoqL4J2jLnTjNqCFhEbgfeUFX0Q+CdJK4C7gQHgX9LpNuoSEe8eZfnngM/Vuz4zM7O8LF88sE9/LoCe7i6WLx7IMara2jHmTtPwIIH0EuY7JM0CFgBPRcRNTY/MzMysDZQ61Z9zzV08N7KH/jYYEdmOMXeaMU2zARAR20imwDAzM+toSxb2c+UtDwCw6oyjR6ldDO0Ycydpxq2ezMzMzKyJnKCZmZmZFYwTNDMzM7OCcYJmZmZmVjBO0MzMzMwKxgmamZmZWcE4QTMzMzMrGCdoZmZmZgXjBM3MzMysYJygmZmZmRWMEzQzMzOzgsklQZN0maRHJd1dZbkkfUbSFkl3SXp9q2M0MzMzy0teLWhfAY6tsfz3gfnpYxnwxRbEZGZmZlYIuSRoEfFD4IkaVU4EvhqJm4BeSTNbE52ZmZlZvoraB60feLDs+da0bD+SlklaL2n99u3bWxKcmZmZ2UQqaoKmjLLIqhgRF0fEoohY1NfXN8FhmZmZmU28qXkHUMVWYE7Z89nAtpxiMTMz22vNhkFWrtvMtqFhZvX2sHzxAEsWZl7kaUuTffvaRVETtLXA2ZKuAn4LeCoiHso5JjMz63BrNgyyYvVGhneNADA4NMyK1Rtzjqp5am2fk7TWymuajSuBHwMDkrZKOl3SmZLOTKtcC9wHbAG+DPxZHnGamZmVW7lu897kpWR41wgr123OKaLmmuzb105yaUGLiHePsjyAs1oUjpmZWV22DQ1XLZ99cE+Lo2m+WttnrVXUQQJmZmaFM6s3OwmrVt5uJvv2tRMnaGZmZnVavniAnu6ufcp6urtYvnggp4iaa7JvXztxgmZmZlanJQv7ufCkBRzQlZw++3t7uPCkBZOmA/1k3752UtRRnGZmZoW0ZGE/V97yAACrzjg652iab7JvX7twC5qZmZlZwThBMzMzMysYJ2hmZmZmBeMEzczMzKxgnKCZmZmZFYwTNDMzM7OCcYJmZmZmVjBO0MzMzMwKJrcETdKxkjZL2iLpYxnLT5O0XdId6eODecRpZmZm1mq53ElAUhfweeBtwFbgVklrI+KeiqqrIuLslgdoZmZmlqO8bvV0JLAlIu4DkHQVcCJQmaCZmZnlYs2GQVau28y2oWFm9fawfPFAx96T0vui9fK6xNkPPFj2fGtaVumdku6SdI2kOVkrkrRM0npJ67dv3z4RsZqZWYdZs2GQFas3Mjg0TACDQ8OsWL2RNRsG8w6t5bwv8pFXgqaMsqh4/h1gXkS8BvgecHnWiiLi4ohYFBGL+vr6mhymmZl1opXrNjO8a2SfsuFdI6xctzmniPLjfZGPvBK0rUB5i9hsYFt5hYh4PCJ2pk+/DBzRotjMzKzDbRsabqh8MvO+yEdeCdqtwHxJh0k6AFgKrC2vIGlm2dMTgHtbGJ+ZmXWwWb09DZVPZt4X+cglQYuI3cDZwDqSxOvqiNgk6XxJJ6TVPiRpk6Q7gQ8Bp+URq5mZdZ7liwfo6e7ap6ynu4vliwdyiig/3hf5yGsUJxFxLXBtRdknyn5fAaxodVxmZra/ThvFV9q2c665i+dG9tDfAdtcTSfvizzf97klaGZmVjxZJySAFas37u0oXhrFt/7+J/j+T7ZP2qRtycJ+rrzlAQBWnXF0ztHka7LviyK+752gmZl1qMqT0u/+Zh/fvG1wvxPSgd1TMkfxXXHTA3uH35fqlnRSa5u1j0YSsbzf907QzMw6UGluq/KTUvmJp2R418h+J6mSrLrnrd3Ezt179jvZAU7SLFdZ7/laiVje73snaGZmk1xWq0HW3FaVJ56xGBretV9Z+ZxZRW1Z67Q+ds1S5P1WGduzz+1uKBFrRLX3/XlrN415/zhBMzObxKq1GjRyUurt6d6ndQCS2cYbSegqX7dILWvV9pHVVmu/FfGYNqoZ7/uh4V17k7dG909e86CZmVmTrdkwyDEX3cBhH/tnjrnohr0tCFmtBl3KuqHL/rd56enu4rwTXsWFJy3ggK7klNHf28N7jjo0c+qFg1/QnbneLqmws9F7pvyxKfJ+y4qtmt6e7sz3cjPe95Ua2T9uQTMzmwQabSkbiaCnu2uf5T3dXbzziH6uvnVr5nQKlaP4Fs09ZL+pF4D9XrfydcoNDg1zzEU35HqJrNZM+bMP9mSs1RTlDgNZl1nrjaGUiEH1aUTG+r6vpt7YnKCZmbWhevvXdEmMxP4XZUonlqyT0n8+sgMYfTqFWlMvVK535brNmZeZxPOXn/K6RDartyczNs+UX1sR9lu1Lya9L+jmyWf37xfW29PNs8+N1JWIVVPv+/7Z53ZnxlDv/nGCZmbWZhrpX1Otpax0YpqIua2qrbeyhSGrP08eAwqWLx7IbPVbvnhg73bY/mrtt4lS7xeTaVOnZL7vzzvhVRM2n1vl+77yc1qKod794wTNrEmqza9T1BFO1h7qHYFZTa2WslbKmo2+WmI5kQMKao06zNpHTtCqq7bfgAm5bN3IF5Onhnfx96e8LtdjOt47MDhB6wDV/iE1Wt6J6k26YP+JDpd/404Q7BqJvWWeyNNqqXfi2HqTs4luKWtUZRzHXHRD5km22oCC8UxZAKOPOizCPmo3o7UajTW5Hu8Xk1m9PYU4puOJwQnaJFPvP/j19z/RcPlkvqVLvbNLV0u6siY63LVn/34/tSY0nOz7uJpGWh7b8UvFeJL8ahPHVutXVqt/TRFVu0RW7SRcbcoCyN6f9ZzgS5dUi7yf2sloIzvH+lkYyxeTdpdbgibpWOD/AF3AJRFxUcXyacBXgSOAx4FTIuIXrY6znTQyM/iVNz+43z/4WuXteEuX8ZwYG026GplTqtqEhu24jxsx3iS4WV8qxpvMNZIkNrJ9We+3avMtVetXNpH9ayZCtUtA1QYUVKr2hafaPq72OW31qMPJrNq+zLps3chnYTJ9MalXLgmapC7g88DbgK3ArZLWRsQ9ZdVOB56MiFdIWgp8Cjil1no3Dj7FMRfdMGkOTi3jnRk8601eq7yRW1u0uiVoIk6MzZpduhGT6bY547nHY7UkuBlfKsabzDXaIj1RSX6tfmXt1meq3gEF1WR94am2j6ud4D1as3mqjezMumzd6GdhsnwxqZeiygl5Ql9UOho4LyIWp89XAETEhWV11qV1fixpKvAw0Bc1An759IPigte+kSlTRN/0aQw9u4udu0eYNrWLOYf0MGP6tIndsBZ5bMdO7nvsl+wpe3NPmaJ9no9GElm7slr5eEyZIn59xgsBePCJ4X2OSVbZjOnTeGzHzrrKe1/QzfYdO/fbF1Mkdo/saep21DK1awp7IvaJQ+lEoOX7s1mxlfZJ1j5qpazjBOz3/iyKet/fpf8hle+t8a63GUqfpxnTp3HPQ08DcPjMF+9dnlXWaPlE1W1kHY/t2MnPtv+SiGDa1C5GIpryma78X9nM/Vnkfd+q2Jpxfqqm9D+m/H1R+r9XxH1RXvZw3xzef+Vn99smSbdFxKKs7c0rQTsZODYiPpg+fy/wWxFxdlmdu9M6W9PnP0vrPFaxrmXAMoBff+GLj7jgtW/MfM3yJCHr4EJxPyz1/qMa78mn0ZNSIxpJYFodWzNibiQBheYkMFknmr7p03j0mZ0teX9X+0c8UclxK79UTOR6q8l6v032L5uNaMb7rShfbCa7rC9upef1qPZZKCXS7WjaK3+Tl5177n7ltRK0vPqgZd1jpPI/YT11iIiLgYsBps2cHx/97T+r+qJZ99Xq6e7iwpMWsGRhP+f844+BfZtIs8qqlU9E3ax5VGqpNjN45eWbN1XpM5NVXnlZp7TeA7unZE7CN17VLkNUKx+vWu+LLurr1HpEernxiIz1Z5U9OM59XG1fVM4rVdoOyB7qnfU+XLNhMHPY/Fj7CdWSte+7p2ifS9Gl7XjnEf2Z+yirvNr98ibqPVSrb0y921ft/famAl/KbqW57P+5Kb03/7Lif2StfXzEwv7Mz6Q1z1z2/7/34IZBzmvgOGV9Fo7osM9CXgnaVmBO2fPZwLYqdbamlzgPAp4Yz4tW65xdGl2y4YEhnhvZk1s/tjUbBveLYSzzHdXbAXrJwv7MZVnli+YeMmo/L2j8RrJZGu0fl6WRE2PpNh/V9lu1fTQe49nHtUa51erH9lzaylDeb6ry/VZ6vfK6jXa4riYreay277PKlizsz9xHWeXVEt5mJHNZ21EteWx0+6DYfQvzVu1/FjS2j631Svven4X65XWJcyrwU+AtwCBwK/BHEbGprM5ZwIKIODMdJHBSRLyr1nqnzZwfM0/9hzElCY3ek+6UjJaHrLJ661abcbiRlrNSS2ArTURrW6MtaLVajSbDP+1qA0LG23qVtd+acZyqtUhmteS2evBIPR3/y+OtVj4Ro0PNrPMUrg8agKTjgH8gmWbjsoj4O0nnA+sjYq2kA4GvAQtJWs6WRsR9tdY5bZgIVxkAAAWySURBVOb8WPThf8ztclG9CVrWJaRqJ9xaJ8EXTptayJNBPaMqoTmXsvI48RdBVkLfjNbLRmV9samWHBf5mLTj/Gpm1v4KmaBNhEWLFsX69euB+pOERlqpoHbrQFZrW2UyVi15rBVDtZNgu50oJvuEpK2Wd1/BRi+pm5nZvjoyQaumlZeLmtHfxSdBa0S9X0yqvQ8b7dTu96GZ2dg5QRvFRF0uanTE2GRpKbPiqbe1bTL13TMzKzonaHXw5SLrNL5cbGaWr45J0CRtB+5v1vqm9Lz4kK7ph/Sra+oBMbL7uZEdTwwCTH1x31ykKaOuIMiezW2/erFn99Pb798z/PRo04jMAB4bpY4Vk49de/Pxa18+du1tsh+/uRHRl7VgUiVok52k9dUybSs2H7v25uPXvnzs2lsnH7/RW4HMzMzMrKWcoJmZmZkVjBO09nJx3gHYmPnYtTcfv/blY9feOvb4uQ+amZmZWcG4Bc3MzMysYJygtQlJfyEpJM1In0vSZyRtkXSXpNfnHaPtT9JKST9Jj9G3JPWWLVuRHr/NkhbnGadlk3Rseny2SPpY3vFYbZLmSPq+pHslbZL04bT8EEnXS/rP9OfBecdq2SR1Sdog6bvp88Mk3Zweu1WSDsg7xlZxgtYGJM0B3gY8UFb8+8D89LEM+GIOodnorgdeHRGvAX4KrACQdDiwFHgVcCzwBUlduUVp+0mPx+dJPmuHA+9Oj5sV127gzyPilcBRwFnpMfsY8G8RMR/4t/S5FdOHgXvLnn8K+Pv02D0JnJ5LVDlwgtYe/h44h33vPnUi8NVI3AT0SpqZS3RWVUT8a0TsTp/eBMxOfz8RuCoidkbEz4EtwJF5xGhVHQlsiYj7IuI54CqS42YFFREPRcTt6e/PkJzo+0mO2+VptcuBJflEaLVImg28HbgkfS7gzcA1aZWOOnZO0ApO0gnAYETcWbGoH3iw7PnWtMyK6wPAdenvPn7F52PUxiTNAxYCNwMvjYiHIEnigF/LLzKr4R9IGiP2pM9fAgyVfcntqM/g1LwDMJD0PeBlGYs+DpwL/F7Wn2WUeUhuDmodv4j4dlrn4ySXX64o/VlGfR+/YvExalOSpgPfBD4SEU8nDTFWZJLeATwaEbdJelOpOKNqx3wGnaAVQES8Natc0gLgMODO9B/MbOB2SUeSfJOYU1Z9NrBtgkO1DNWOX4mkU4F3AG+J5+e18fErPh+jNiSpmyQ5uyIiVqfFj0iaGREPpV1BHs0vQqviGOAESccBBwIvJmlR65U0NW1F66jPoC9xFlhEbIyIX4uIeRExj+SE8fqIeBhYC7wvHc15FPBUqQnfikPSscBHgRMi4tmyRWuBpZKmSTqMZLDHLXnEaFXdCsxPR5EdQDKoY23OMVkNaZ+lS4F7I+LTZYvWAqemv58KfLvVsVltEbEiIman57qlwA0R8R7g+8DJabWOOnZuQWtf1wLHkXQufxZ4f77hWBWfA6YB16etoDdFxJkRsUnS1cA9JJc+z4qIkRzjtAoRsVvS2cA6oAu4LCI25RyW1XYM8F5go6Q70rJzgYuAqyWdTjIa/g9zis8a91HgKkl/C2wgScA7gu8kYGZmZlYwvsRpZmZmVjBO0MzMzMwKxgmamZmZWcE4QTMzMzMrGCdoZmZmZgXjBM3MzMysYJygmZmZmRWMEzQzsyok/ZOkd0jqlXSdpD/IOyYz6wxO0MzMqlsADJHcXuZvIuJbOcdjZh3CdxIwM8sgaQrwDPA48PmI+FTOIZlZB3ELmplZtvnANuA04ExJ3fmGY2adxAmamVm2BcD1EXEDcDfwvpzjMbMO4gTNzCzbApLEDOACYIWkqTnGY2YdxH3QzMzMzArGLWhmZmZmBeMEzczMzKxgnKCZmZmZFYwTNDMzM7OCcYJmZmZmVjBO0MzMzMwKxgmamZmZWcE4QTMzMzMrmP8PScmox/NRFRgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAADgCAYAAABsKDD3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5xU9X3v8ddnZmdh+bkgqOyCgIo/UMTVwUS9zY2mCcabKrWaaKImTRqT3PxomoRE2t5cb35pQlvbNMbGmzRNtPFnEGk0JYna9iaNuosgCIggGmQHkF+L/Fj25+f+cc5sh92Z3Rl2ZmfO7vv5ePBgzznfOedz5pyZ85nvj3PM3RERERGRaImVOwARERERKZySOBEREZEIUhInIiIiEkFK4kREREQiSEmciIiISAQpiRMRERGJICVxIsOMmf2emW0qdxzZmNnbzWz7EG3rn8zsa0OxrXyY2c/N7INl2O6HzOzXQ71dESk9JXEiFcLMXjOzVjM7lPHvO3m8zs3s9PS0u/8/dz+zRDFWVGJ0vMqR2Lj7u939R0O5TREZ3qrKHYCIHOMP3P1X5Q5CwMzi7t5V7jiiwMyq3L2z3HGIjDSqiROJADM73cz+3cwOmNkeM3swnP8fYZEXwpq79/Vusgxr+Bab2VozO2xmPzCzk8LmvYNm9iszm5RR/mEz2xlu6z/M7Jxw/i3AB4Avhtv6l3B+nZn91Mx2m9mrZvaZjHXVhLV3+81sA7BggP38OzN73czeNLNVZvZ7GctuM7OHzOzHYdzrzSyZsbzBzJ4Plz0IjM6xjbOBfwAuDvejJZz/T2Z2t5k9YWaHgcvM7N/M7E8yXntMDZ6ZnWVmvzSzfWa2ycze28++9awrXM9vzOxOM2sxs61mdkk4/3UzeyOz6TWM7R/CbR0Mz4WZ4bJZYW1sVbZt9YrBwm2+ER7ftWZ2brhslJn9lZltM7Nd4fZqcuxLZvz7gNvC+R82s43h8V6ZEWN/2825b+HyS8ysMXxdo5ld0ms/vxrGctDMfmFmU8Jlo83sPjPbG77HjWZ2UrhsogWfgx1m1mxmXzOzeK5jJ1KplMSJRMNXgV8Ak4DpwN8DuPvbwuXz3X2cuz+Y4/V/BLwTOAP4A+DnwJ8DUwi+Bz6TUfbnwBzgROB54J/Dbd0T/v2tcFt/YGYx4F+AF4B64B3AZ81sYbiu/w2cFv5bCAzUJ6wROB+YDPwEeNjMMpOxq4AHgFpgBfAdADOrBpYD94avfTjc5z7cfSPwceC34X7UZix+P/B1YDzQb3OrmY0FfhnGeSJwA/BdC5PePLwFWAucEK7jAYIk93TgRuA7ZjYuo/wHCM6DKcAawuNSoHcBbyM4D2qB9wF7w2XfDOefH8ZQD3x5gPi3Euz7181sEcE5dQ0wFfh/wP15bDfnvpnZZOBx4NsE79PfAI+b2QkZr30/8MdhHNXAF8L5HwQmAjPC134caA2X/QjoDPezIYyvT9IrUumUxIlUluVhrUH630fD+R3ATKDO3Y+6e6H9uf7e3Xe5ezPBxfVZd1/t7m3AowQXMgDc/R/d/WC47DZgvplNzLHeBcBUd/+Ku7e7+1bg/wLXh8vfC3zd3fe5++sEF+Oc3P0+d9/r7p3u/tfAKCCzf9+v3f2JsJnzXmB+OP+tQAL4W3fvcPdHCBLCQj3m7r9x9253PzpA2fcAr7n7D8N4nwd+Clyb57ZeDV/bBTxIkGx8xd3b3P0XQDtBkpH2uLv/R3hc/oKgJnFGQXsXnEfjgbMAc/eN7r7DzAz4KPBn4bE6CHyD/zqO2aTc/e/DfW8FPgbcHq6zM3z9+WGtWtbt5rFv/wPY7O73htu5H3iJ4IdI2g/d/eUwhocIktD0vp4AnO7uXe6+yt3fDGvj3g181t0Pu/sbwJ0D7KtIRVISJ1JZFrl7bca//xvO/yJgwHMWNCN+uMD17sr4uzXL9DgI+oGZ2R1m9oqZvQm8FpaZkmO9M4G6zMSToDbmpHB5HfB6Rvnf9RekmX0+bI47EK5rYq9t78z4+wgwOmxGrAOa3d3z3VYOrw9cpMdM4C299v0DwMl5vr73McDdsx6X3rG5+yFgH8F+583dnyKovbwL2GVm95jZBIKaszHAqox9+ddwfi6936uZwN9lvH4fwTlb3892B9q3Ovoex98R1BKm9T4n0u/ZvcBK4AEzS5nZt8wsEcaZAHZkxPo9gpo8kUhREicSAe6+090/6u51BDUe37WMEalF9H7gauD3CRKoWeF8S4fSq/zrBDVKmYnneHe/Mly+g6CGKe2UXBu2oP/blwhq7yaFzZwHMrbdnx1AfVijNOC2suxHrvmHCZKbtMwE7XXg33vt+zh3/0Qe8R6PnvcxbGadDKTCGOknzmO4+7fd/ULgHILmzcXAHoKk8ZyMfZno7uNyrYfs58LHer0fNe7+n/1sd6B9SxEkXZlOAZr7iSu9nx3u/n/cfS5wCUHN6c1hnG3AlIw4J7h7vs3gIhVDSZxIBJjZdWY2PZzcT3ABTY+c3AWcWqRNjSe4wO0lSAq+0Wt57209B7xpZl+yYBBD3MzONbP0AIaHgCVmNimM/9MDbLsT2A1UmdmXgQn9lM/02/C1nzGzKjO7Brion/K7gOlhX7r+rAGuMbMxYdL8kYxlPwPOMLObzCwR/ltgwcCJUrjSzP5bGPNXCZrEX3f33QRJzY3h+/9hgj6IfYTxvSWskToMHAW63L2boBn8TjM7MSxbn9G3MR//QHCs0wNhJprZdf1td6B9A54geI/fHx7X9wFzCd77fpnZZWY2Lxyw8CZB82pX2Iz7C+CvzWyCmcXM7DQz++8F7KtIRVASJ1JZ/sWOvU/co+H8BcCzZnaIoEP/n7r7q+Gy24AfhU1DOUdH5unHBM1VzcAG4Jley38AzA23tTzsz/UHBP2QXiWo0fk+QS0ewP8J1/cqwYXz3n62vZJgUMXL4WuOkmfzpru3E3So/xBBkvs+YFk/L3kKWA/sNLM9/ZS7k6Bv2i6CzvA9gwnCfmPvIuhLlSJo1vsmQT++UvgJwUCRfcCFBE23aR8lqNnaS1DT9Z851jGBIFnbT/Ae7wX+Klz2JWAL8EzYlP4rju2P2C93f5Rg/x8IX/8iQd+zgbabc9/cfS9BDdrnw9d8EXiPu/d3zNJOBh4hSOA2Av8O3Bcuu5lgEMSGMKZHgGn57qtIpbBju5CIiEilMbN/Ara7+1+WO5ZiG877JlJqqokTERERiSAlcSIiIiIRpOZUERERkQhSTZyIiIhIBCmJExEREYmgqoGLDD9TpkzxWbNmlTsMERERkQGtWrVqj7v3eYLKiEziZs2aRVNTU7nDEBERERmQmWV9jKCaU0VEREQiSEmciIiISAQpiRMRERGJICVxIiIiIhGkJE5EREQkgpTEiYiIiESQkjgRERGRCBqR94kTEZHyW766maUrN5FqaaWutobFC88E6DNvUUN9mSMVqUxK4kREZMgtX93MkmXraO3oAqC5pZXFD78ABh1d3jNvybJ1AErkRLJQc6qIiAy5pSs39SRwaR3d3pPApbV2dLF05aahDE0kMlQTJyIiQy7V0ko3R9g1agld1tJneaJ7Fie234ZhpFpayxChSOVTEiciIkOurraGLW+uoT32CjVdbyHuE3uWdViKo/FVdLGXKqZQV1tTxkhFKpeSOBERGXKLF57JLY/eB8AJ7X9GnHEkYgYGh7o3sjP+BdpjmxkfP6lnwIOIHEt94kREZMgtaqjntOk7Sfg04oyjvraGpdfNZ+m18xkbOw08TvWYV7n9mnka1CCSg2riRESkLFJH1nNyzTm8ZdpkHvzYxT3z739uG3v3nMo5M/cogRPph2riRERkyL1x+A22HdjG5MRZWZdPSpxFU6oJd8+6XESUxImISBk0pZoAmJw4O+vyyYmz2de6j1dbXh3KsEQiRUmciIgMucbmRgyjtuqMrMvTyV062RORvpTEiYjIkGva0cRZU84iERubdfmEqlOpjlfT2Nw4xJGJRIeSOBERGVLuTlOqiQX1C3KWiVuC808+n6YdqokTyUVJnIiIDKnmg83sPLST5LRkv+WS05KsSq2i27uHKDKRaKnoJM7MrjCzTWa2xcxuzbL8FDN72sxWm9laM7uyHHGKiEj+0v3cknUDJHF1SQ62H+TlvS8PRVgikVOxSZyZxYG7gHcDc4EbzGxur2J/CTzk7g3A9cB3hzZKEREpVGNzI3GLc/7J5/dbLt3cqsENItlVbBIHXARscfet7t4OPABc3auMAxPCvycCqSGMT0REjkPTjibOPfFcahL9PxP1rClnMSYxRoMbRHKo5CSuHng9Y3p7OC/TbcCNZrYdeAL4dK6VmdktZtZkZk27d+8udqwiIpKHnkENdbkHNaRVxaq4YNoFGtwgkkMlJ3GWZV7vW3ffAPyTu08HrgTuNbOs++Tu97h70t2TU6dOLXKoIiKSj1dbXmVf674B+8OlJaclWb1jNZ3dnSWOTCR6KjmJ2w7MyJieTt/m0o8ADwG4+2+B0cCUIYlOREQKlu+ghrRkXZLWzlY27N5QyrBEIqmSk7hGYI6ZzTazaoKBCyt6ldkGvAPAzM4mSOLUVioiUqEamxupjlcz76R5eZXX4AaR3Co2iXP3TuBTwEpgI8Eo1PVm9hUzuyos9nngo2b2AnA/8CHX05JFRCpW044m5p80n+p4dV7lT598OhNGTdDgBpEsqsodQH/c/QmCAQuZ876c8fcG4NKhjktERArX7d2sSq3ipvNuyvs1MYuRrEtqcINIFhVbEyciIsPLy3tf5mD7wbz7w6UlpyV5YecLtHW2lSgykWhSEiciIkOi0EENacm6JB3dHax7Y10pwhKJLCVxIiIyJBqbGxmTGMPZU88u6HUa3CCSnZI4EREZEk07mmg4uYGqWGHdsWdOnMkJNScoiRPpRUmciIiU1PLVzVxy+y/57bZVbE2dzPLVzQW93sw4Zfw87l/zNLNvfZxL73iq4HWIDEcVPTpVRESibfnqZpYsW8eBzi346DbaW2exZFnQt21RQ+8nKeZex7ad0zhi/0EXR2luoeB1iAxHqokTEZGSWbpyE60dXbTFNgMwqnsOrR1dLF25qaB1xDpOA+umPfYqQMHrEBmOlMSJiEjJpFpaAWiPbcZ8DFVed8z8fNdR3T0nWI+93GfdIiOVkjgRESmZutoaANpjW6juPh0LLzvp+fmuo4oTiPtk2mNb+qxbZKRSEiciIiWzeOGZjEp00W6vMiqsTatJxFm88MyC1lGTiFPdPaenWbbQdYgMR0riRESkZBY11PPRy6vAOqnunkN9bQ23XzOvoAEJixrquf2aedT4HDqtmZMmdhe8DpHhSEmciIiU1Ljx2wC4qD7Jb269/LiSr0UN9ZwztQHM+dYNE5TAiaAkTkRESqyxuZFqm8jY+LRBrWdS4qxgfanGYoQlEnlK4kREpKSadjQxOXEWZjao9YyOTWJM7GQ9uUEkpCRORERK5kjHEda/sZ5JicKel5rL5MRZqokTCSmJExGRklmzcw1d3sXksCl0sCYlzmbr/q3sa91XlPWJRJmSOBERKZl00+fkotXEBetZlVpVlPWJRJmSOBERKZnGVCPTxk2jJj61KOublDizZ70iI11FJ3FmdoWZbTKzLWZ2a44y7zWzDWa23sx+MtQxiohIbk2pJpJ1yaKtrzo2njmT52hwgwgVnMSZWRy4C3g3MBe4wczm9iozB1gCXOru5wCfHfJARUQkqzfb3mTTnk0sqFtQ1PUm65KqiROhgpM44CJgi7tvdfd24AHg6l5lPgrc5e77Adz9jSGOUUREcnh+x/M4XtSaOIAFdQvY/uZ2dh7aWdT1ikRNJSdx9cDrGdPbw3mZzgDOMLPfmNkzZnbFkEUnIiL9Sjd5FjuJS69PgxtkpKvkJC7bXSG913QVMAd4O3AD8H0zq826MrNbzKzJzJp2795d1EBFRKSvxlQjMyfOZOrY4gxqSGuY1kDMYmpSlRGvkpO47cCMjOnpQCpLmcfcvcPdXwU2ESR1fbj7Pe6edPfk1KnF/UIREZG+ij2oIW1c9TjOnnK2BjfIiFfJSVwjMMfMZptZNXA9sKJXmeXAZQBmNoWgeXXrkEYpIiJ97Gvdx9b9W4s+qCEtWZekKdWEe+8GGpGRo2KTOHfvBD4FrAQ2Ag+5+3oz+4qZXRUWWwnsNbMNwNPAYnffW56IRUQkrVT94dIW1C1g1+FdbH9ze0nWLxIFVeUOoD/u/gTwRK95X87424HPhf9ERKRCpJO4C+suLMn608lhU6qJGRNnDFBaZHgqqCbOzCbn8S/rwAIRERk5GlONzJk8h9rRpbkkzD95PlWxKg1ukBGt0Jq4VPgv28jRtDhwynFHJCIikdeUauJtM99WsvWPrhrNvBPnaXCDjGiFJnEb3b2hvwJmtnoQ8YiISMTtPLST7W9uJzmtNP3h0pJ1SR7Z8Ajujll/dQsiw1OhAxsuLlIZEREZpko9qCFtQd0C9h/dz9b9uimBjEwFJXHufhTAzL7We1n4rNOeMiIiMjI1pZqIWYyGaf023Axa5uAGkZHoeG8xUm9mN6QnzOxE4FfFCUlERKJo+epmLr3jKb711BOM4hR+tf5ASbf3SmoSRoJbHniIS+94iuWrm0u6PZFKc7y3GPkYsNLMXiF4FNYPgS8VLSoREYmU5aubWbJsHUc6OmkbvZmajiRLlq0DYFFD78deF2d7/2v5SyRip9IWe5nmltaSbk+kEhV6i5Efm9lnCfq9fRK4B7gbWOTuj5cgPhERiYClKzfR2tFFl+2h21oY1X06rR1dLF25qaTbG9V9Ou2xV3C6S7o9kUpUaHPqj8LXfBj4CTAL2A/caGbXFjc0ERGJilRLKwDtthmA6u45x8wv1faqu+fg1kqnNZd0eyKVqKDmVHd/EngyPW1mVcBcYD7wVuCRokYnIiKRUFdbQ3NLK22xzeBxqn12z/xSbi+dLLbFNpPomlGy7YlUokE9O9XdO919rbvf6+5fKFZQIiISLYsXnklNIk57bAvVPgujmppEnMULzyzp9hI+HfPRtMc2l3R7IpWo0D5xzxejjIiIDC+LGur5xh+eS3tsM9Xdc6ivreH2a+aVbJDBooZ6br9mHqPiCaq7T4PEKyXdnkglKnR06tlmtraf5QZMHEQ8IiISUefNaqPbDnHm5Pn85s8uL/n2FjXUc/9z27A3z2Vb+2O8Z/5JJd+mSCUpNIk7K48yXccTiIiIRFv6YfSTE/lcKopncuJsXj7yIOvfWM/8k+cP6bZFyqnQgQ2/AzCzy4EPAC3Ai8Ba4EV3byt6hCIiEglNqSZiVDOx6rQh3e6kxNk921cSJyPJ8Q5suA/4GfAMcCrwZWB9sYISEZHoaUw1UpuYQ8yO9z7yx2dcvJ6Joyb21ASKjBTH+0nb4u6Phn8/XKxgREQkmrq6u3h+x/OcnLhiyLdtFiNZl9QzVGXEOd4nNvzWzD5fophERCRiXt77MofaDw15f7i0ZF2StbvW0tapXj0ychzvExtOBm4ys9+Z2Qoz+6qZXVfs4MzsCjPbZGZbzOzWfspda2ZuZslixyAiIgNLN2VOqipPEregbgEd3R2s3dXfDRREhpdiPrHhIorYtGpmceAu4J3AdqDRzFa4+4Ze5cYDnwGeLda2RUSkME2pJsYmxjK+amZZtp+sS/bEsaB+QVliEBlqxXxiw+JiBRW6iKDv3VZ3bwceAK7OUu6rwLeAo0XevoiI5Kkx1cgF0y4gZvGybP+UiacwZcwUDW6QEWVQSVyJ1QOvZ0xvD+f1MLMGYIa7/2woAxMRkf/S0dXBmp1remrDysHMWFC3QIMbZESp5CTOsszznoVmMeBOIK8BFmZ2i5k1mVnT7t27ixSiiIhs2L2Bo51HWVBX3mbMZF2S9bvXc6TjSFnjEBkqlZzEbQdmZExPB1IZ0+OBc4F/M7PXgLcCK3INbnD3e9w96e7JqVOnlihkEZGRJ92EWc6aOAgGN3R7N6t3rC5rHCJDpShJnJlNM7NRxVhXhkZgjpnNNrNq4HpgRXqhux9w9ynuPsvdZxHcePgqd1dduojIEGpKNTFx1EROn3x6WePIHNwgMhIUqybuXuAlM/urIq0Pd+8EPgWsBDYCD7n7ejP7ipldVaztiIjI4DSmGknWJTHL1gtm6EwbP4368fUa3CAjRkG3GDGzyTkWvReIE9Sc1bp7y6AjA9z9CeCJXvO+nKPs24uxTRERyd/RzqOs27WOz138uXKHAqAnN8iIUuhjt1Lhv94/t9IDDowgmTtlkHGJiEgErN21lo7ujrIPakhL1iV5bNNjHDh6gImjJ5Y7HJGSKjSJ2+juDf0VMDP1KBURGSHStV7lHtSQlk4mn9/xPJfNvqzM0YiUVqF94i4uUhkRERkGmlJNTBkzhVMmVkYDzIV1FwIa3CAjQ0FJnLsfBTCzr/VeFj4mq6eMiIgMf42pRhbULSj7oIa0KWOmMLt2tgY3yIhwvKNT683shvSEmZ0I/Ko4IYmISBQcbj/Mht0bKqYpNU2DG2SkON4k7mPALWZ2kZktAJ4CinZ7ERERqWzLVzfz1r/6Pt3ezcP/Wc3y1c3lDqlHjc/h1ZZXOeXWn3DpHU9VVGwixVToLUZ+DDwPrAY+CfwE6AQWufuW4ocnIiKVZvnqZpYsW8eu7vVQDQcPnsKSZevKHRYQxPb0ixMgDm2xLTS3TOyJbVFD/QCvFomWQmvifhS+5sMECdwsYD9wo5ldW9zQRESkEi1duYnWji7aY5uJ+2SqOIHWji6WrtxU7tCCGNpnA9AeC+oWKiU2kWIrqCbO3Z8EnkxPm1kVMBeYT/Ds0keKGp2IiFScVEsrAG2xzVR3zzlm/vRJNeUKqyeGGGOp6p5OW2zzMfNFhptBPXbL3Tvdfa273+vuXyhWUCIiUrnqamvo5jCdseZjkri62vImcJkxjOo+nfaMJK4SYhMptoKSODN7vhhlREQkuhYvPBOqtwIwKkziahLxYH6ZLV54JjWJONXdc+iyvXSyr2JiEym2Qp/YcLaZre1nuQF6zomIyDC2qKGeRzcf5Mcbobr7dOpra1i88EwWNdRz/3Pbyh4bwKd/+iL7gfHjf8ftV16uQQ0yLBWaxJ2VR5mu4wlERESio9U2MzY+jUtmz+bBj1XWg3oWNdRz7zMX0PxGjOsuaVcCJ8NWoQMbfgdgZpcDHwBagBeBtcCL7t5W9AhFRKTiNKYamZQ4u9xh5FQVq2FC1Wzd9FeGtX77xJnZXDO7L8ui+4CfAc8ApwJfBtYXPzwREak0e47s4bWW15hclU/jTPlMSpxFU6oJdy93KCIlMVBN3JNkf6D9Fnd/NPz74eKGJCIilWxVahUAkyu4Jg5gcuIsXnvzcbYd2MbM2pnlDkek6AYanfou4OvpCTP7sZl9FvitmX2+pJGJiEhFSj9cvjZR2SM+00mmmlRluOo3iXP3de7+gYxZ6Sc2nAzcZGa/M7MVZvZVM7uulIGKiEhlaEo1ccYJZ1AdG1fuUPo1sep0ErFET9IpMtwU84kNF6GmVRGRYa8x1chlsy6jY3e5I+lf3Ko576TzVBMnw1Yxn9iwuFhBpZnZFWa2ycy2mNmtWZZ/zsw2mNlaM3vSzNTpQUSkhFIHU6QOpkjWJcsdSl6SdUkNbpBha1BJXCmZWRy4C3g3QW3fDWY2t1ex1UDS3c8jeG7rt4Y2ShGRkSU9qGFB3YIyR5KfZF2SA20H2LJvS7lDESm6ik3iCJpnt7j7VndvBx4Ars4s4O5Pu/uRcPIZYPoQxygiMqI0phqJWYzzTz6/3KHkJZ1sqklVhqNKTuLqgdczpreH83L5CPDzXAvN7BYzazKzpt27K7wjh4hIhWpKNXHO1HMYWz223KHkZe7UuYyuGq3BDTIsVXISZ1nmZe3UYGY3Aklgaa6Vufs97p509+TUqVOLFKKIyMjh7jSlmiLTHw4gEU9w/snnqyZOhqVKTuK2AzMypqcDqd6FzOz3gb8ArtJjv0RESmfbgW3sPrI7UkkcBE2qz+94nq5uPdpbhpdKTuIagTlmNtvMqoHrgRWZBcysAfgeQQL3RhliFBEZMdK1WVEZ1JCWrEtyuOMwL+15qdyhiBRVxSZx7t4JfApYCWwEHnL39Wb2FTO7Kiy2FBgHPGxma8xsRY7ViYjIIDWmGknEEpx30nnlDqUgGtwgw1VBN/sdau7+BPBEr3lfzvj794c8KBGREaop1cR5J53HqKpR5Q6lIGeccAbjqsfRmGrkg+d/sNzhiBRNxdbEiYhI5YjioIa0eCzOBdMuUE2cDDtK4kREZEBb9m3hQNuByPWHS1tQt4A1O9fQ0dVR7lBEikZJnIiIDChdixXFmjgI4m7rauPFN14sdygiRaMkTkREBtSUamJ01WjmTu399MNo0OAGGY6UxImISE7LVzdz6R1PcddvVpLoOpXH10bzbk6nTjqVsYmJ/K8nHmP2rY9z6R1PsXx1c7nDEhkUJXEiIpLV8tXNLFm2ju0th2iPvQIdp7Fk2bpIJj+PrUnhbaeyv/MlHGhuaY3svoikKYkTEZGslq7cRGtHFx3WjNtRqrvn0NrRxdKVm8odWsGWrtxEVdfptNtrOO0Akd0XkTQlcSIiklWqpRWA9liQ6IzqnnPM/ChJtbQG8VsXbbFXjpkvElUVfbNfEZFSW766maUrN5FqaaWutobFC88E6DNvUUN9mSMdenW1NTS3tHIo/jTx7hOp8vqe+VFTV1vDtpZ54AkOx59mdPfZPfNHmmzn/KKG+pzzpXIpiRORESFXsrZk2TpaO4IHoze3tLL44RfAoKPLe+YtWbauZz0j6SK3eOGZfG7ZE7TF11LbcTNGjJpEvOe9i5LFC89kybJ2xnb9HofjTzOp40OMTYyP5L4Uovd5f9lZU/npquZjzvkly9bR9Lt9WeenjaTzPkqUxInIsJfuoN/7AjU6EeuZl9bR7X1e39rRxW0r1tPW2Z31IjdcL2iLGur5wYvP8uorVYzrfBf1Eb6Ap2P+9E/fw2F7isT4/+T2K78QyX3JV7bz/p+f2UbvM7y1o4v7n32dLvc+80fieR8lSuJEZNjpXftwpL2zT7LW2tHVZ15/Wlr73ljm994AABABSURBVOk/s2P8cKypONJxhF+nfsqMmrdz8bRTefBjF5c7pEFZ1FDPT559C7/cO4fxJzzN1ef/TblDKppsNc3pgSmZ+v5ECfRO4NJynfe3rVg/LM/5qNHABhEZVtK1D80trT23kth/pHSPWkrXTGRub7jcuuLBFx+k5WgLp4/5w3KHUjRmxmlj/pAXdr3AM9ufKXc4RZHtnE9P5ytuVtA2W1o7huU5HzWqiRORyMq39iGX2prEMU1FAImYHdMnDqAmEWd0IpY1GYybZa3lW7pyU+RrJu5uupu5U+cyJXZ+uUMpqlNGv4uXj97N3U13c/GMaNcuAlnP+daOLuJmWWvYjGNr5GoScf7owvpj+sSl5+c673sb7rXSlUpJnIhEUq5+bvkmcDWJOLdddQ4AX3xkLe1d3T19vnLN673+mkQ85/aaW1q59I6nInsxW5VaRWOqkW9f8W1+vaawWppKl4iN4abzbuIHq3/AnQvv5IQxJ5Q7pIL0/vGSq8aty73POZpO2B5q3H7M+b2ooZ7kzMl5nfe59P4Mqv9c6SmJE5GKl2+NW3+1D7U1CY60d/W5cAHc/9w2gGP6fGWbB32Tu6UrN2W9iBr0zI/ixezuprsZkxjDzfNv5tdrNpQ7nKL7ePLjfLfpu/xwzQ/5wiVfKHc4ecv246V3zVpa+hztfc4uaqhn865DwLHn96KG+rzO+yPtnQXVSqv/XOmMyCRuXfMBLr3jKZ1IeSrkPlq655YUW6E1brlqH2676pycF6h85brI9Y4n20U1Ss1NLUdbuP/F+3n/ue9n4uiJ5Q6nJOadNI9LZ1zK91Z9j89d/DliVnldxAsZrJCtiTR9bhX7vO/9mUxvL9dnsqW1o2eARBR/0FSyEZnEwbH3xnn6pd0V/YVabIXc6BHyv49WtvsMjeR7bin5PT75jizNVePWX+1D+kJUTOljkrm9XM1buZqbKu176N4X7uVIxxE+seATZYthKHwi+QlufPRGntz6JO887Z1ljSXf+7nlSpQcqI7HstY0F1u2c76/WuneovSDptiKfaPlik7izOwK4O+AOPB9d7+j1/JRwI+BC4G9wPvc/bV819/a0XXMPXOGY4Ix2Bs9FnIfrWz3GRoJ99wa7slvoUn/YMpme99yyVXjVqzah0L03t6ldzyVNfZczU25vocKOa7FOk7f+teXaDq6lHHxs9i28yQumHYcb0hEXDv3Wv7n45/h+p98hQmH24t6fhd67PK9n1t/P16mTwqePlGOcz6tGP3nIPt3Qym+h4byu/N4brQ8UHwVm8SZWRy4C3gnsB1oNLMV7p7ZOeMjwH53P93Mrge+CbyvkO1k+5BUSoIx2C9l6HtBLPRGj4XcRyvXfYayieo9t/L9UJYj+YXif/FB33Oo0MSzkLLZ3rdc+qtxK7fg6QD5NzcV0vSaa16xjtP+rjV0jHqdCW1/esy5NRz9fN0eqo5ezj5bRg17aG6ZUpTzu7/a1WLczy3Xj5dS1DQXIlsN3fH0n+v9HVeM75ZiHKehvv7mM8K9YpM44CJgi7tvBTCzB4Crgcwk7mrgtvDvR4DvmJm5F5BNZFHozQ2L8UtgMBl7IRfEQm/0WIhcvxALUcnNTYX8Wh7q5LdUX3zZzqFCE89CyhYysrQcNW65HGw7yN724KvpueY4dSfCH1/WyXee2kJHVzdTx4/mxreewn3PbGP3waN5rXPrAfjssnW0dXaBBdOffiQ4Tp1d3jPvs8vWUV1lHOzsDDpGhdo64YfPbaTb/dj5zn91ospS9s3Eo8R8LGO6fo/WruBCkq7hGW6WrtzE6PaFMPoR3kz8lLGdbz/u9633/H98dkPw3ZBxnB7bOJWnN+4+5phmHuN8pM+l3udW3YnNPedgR1eSRDwxyHfn+BSj/1xvxfhuyZUoZasF7y/hG8rrbyqPpmkbZL5TMmZ2LXCFu/9JOH0T8BZ3/1RGmRfDMtvD6VfCMnv6W/dp4yb6N+b/t0HHGIsZU8eNYvehNrozThwLb5qY+d7mKptrfi5mRimOWa71VsVjdLsf9/7lKhszo7OrO+84eovFjFOnjGXKuFF57F3h9hxq4/V9rbR1djGqKs6MyTU908VWqmNaKdsrhqp4jK5ux917jkf62G/Y8SYAc6dNOOY12eaXquxzo3fx3um/OM69qzxXvjaLm1+e2zM9oSZICIr9vhVjHYPZ3jNb9wJwR0Mja6buZrhYtfU6Tuge3TNd7uO051Abr+w+fMznt1Tfp8WQ6ztyKK6/WyfW873zrgaC1obf3Hp5uswqd0/2eW2lfpmb2XXAwl5J3EXu/umMMuvDMplJ3EXuvjfL+m4BbgE4deyEC//6wv9O7ZhE1qQqV4KRI868D+pQnxj56i/BPHXKWIA+CU22eVPGjcqa/OR6/dY9h/tsL59ENi29rt5fDlPGjRrUl0sxYhvq5HeoFeO8zybb+1bqhL0YWmJtrB7d72/HHgdaO9hzqI2Orm4S8RhjR1VxoLXjmPeoWN8Jx3OczI2z909mVHccCD5nDafUDjqWSrR6WwttnV0cTLSzZWJLz/xSnd+FrMfMmFiT4HBbZ8+5MmXcKCbWDFzDdsmRkxlFfNAxldKeQ21Zv2eH+vpbbr2vAekkriYR5/Zr5vW0OOVK4iq5OXU7MCNjejqQylFmu5lVAROBfdlW5u73APcAJJNJX/TvK4DcndL/soCbhg6lQposs92NPn2jx95Nk2/P0dR7YXgCXZhl/dnmzSyg7Os5+oYU8qiYQm5kuXx18zF9NXo3Wadf39+d+fO9+/nt18wjTn79mXK999nKQt9zs5C7qve3H9n0dw71fu9yPemgkLK53rcLK6CfW39mAvMH8fpc/aTy/SyU6jilj8nMCn//j9fq1c3cluXzNNj3Ldd92wYaUV2pfYFLYSbZrwHQ9zuuGN8txThOQ3X9LaR/byXXxFUBLwPvAJqBRuD97r4+o8wngXnu/vFwYMM17v7egdadTCa9qamp3zLZbnNQyIU9m0JPjHwfjdLfBREqd5BANtn6TxT6Qcv3fcu13v7kShqHsr9ePqNhoXhJFZRmVFiushLI9lkox3Ea7sekFKOvc/1A7C9RHu7vcyFK9d0y2ONUzutv5JpTAczsSuBvCW4x8o/u/nUz+wrQ5O4rzGw0cC/QQFADd316IER/8kniesvVObMYv3T7m5/vqBkYPhfEfD9ohdSUFmPQRaX/WlZSNfzoOEVXse8HJqVR6HEq12cykklcqRxPEgelu1+WPtgDG2xzUyFyVYPr17KIiJSDkrgMx5vESWUpVdNrFJuhRURk+IriwAaRfqWTqsH2ccjVn01Jm4iIVDIlcRJpixrq+yRbyZmTs9ai5ZovIiISRWpOFREREalguZpTY+UIRkREREQGR0mciIiISAQpiRMRERGJICVxIiIiIhGkJE5EREQkgpTEiYiIiESQkjgRERGRCFISJyIiIhJBSuJEREREIkhJnIiIiEgEKYkTERERiSAlcSIiIiIRpCROREREJIKUxImIiIhEUEUmcWY22cx+aWabw/8nZSlzvpn91szWm9laM3tfOWIVERERKYeKTOKAW4En3X0O8GQ43dsR4GZ3Pwe4AvhbM6sdwhhFREREyqZSk7irgR+Ff/8IWNS7gLu/7O6bw79TwBvA1CGLUERERKSMKjWJO8nddwCE/5/YX2EzuwioBl4ZgthEREREyq6qXBs2s18BJ2dZ9BcFrmcacC/wQXfv7qfcLcAt4eQhM9tUyHYiZAqwp9xByHHT8Ys2Hb/o0rGLtuF+/GZmm2nuPtSBDChMsN7u7jvCJO3f3P3MLOUmAP8G3O7uDw9xmBXJzJrcPVnuOOT46PhFm45fdOnYRdtIPX6V2py6Avhg+PcHgcd6FzCzauBR4MdK4ERERGSkqdQk7g7gnWa2GXhnOI2ZJc3s+2GZ9wJvAz5kZmvCf+eXJ1wRERGRoVW2PnH9cfe9wDuyzG8C/iT8+z7gviEOLQruKXcAMig6ftGm4xddOnbRNiKPX0X2iRMRERGR/lVqc6qIiIiI9ENJ3DBjZl8wMzezKeG0mdm3zWxL+HiyC8odo/RlZkvN7KXwGD2a+fQRM1sSHr9NZrawnHFKdmZ2RXh8tphZtifMSAUxsxlm9rSZbQwf3fin4fwBH/kolcHM4ma22sx+Fk7PNrNnw2P3YDj4cdhTEjeMmNkMgoEg2zJmvxuYE/67Bbi7DKHJwH4JnOvu5wEvA0sAzGwucD2Qfrzcd80sXrYopY/weNxF8FmbC9wQHjepXJ3A5939bOCtwCfDY5bPIx+lMvwpsDFj+pvAneGx2w98pCxRDTElccPLncAXgcyOjlcT3IbF3f0ZoDa8955UEHf/hbt3hpPPANPDv68GHnD3Nnd/FdgCXFSOGCWni4At7r7V3duBBwiOm1Qod9/h7s+Hfx8kSAbqyeORj1J+ZjYd+B/A98NpAy4HHgmLjJhjpyRumDCzq4Bmd3+h16J64PWM6e3hPKlcHwZ+Hv6t41f5dIwizMxmAQ3AsxT4yEcpm78lqLBIP6XpBKAl44fwiPkMVuQtRiS7AR5V9ufAu7K9LMs8DUkug/6On7s/Fpb5C4Kmnn9OvyxLeR2/yqJjFFFmNg74KfBZd38zqNCRSmZm7wHecPdVZvb29OwsRUfEZ1BJXIS4++9nm29m84DZwAvhl9B04Hkzu4jgF8mMjOLTgVSJQ5Usch2/NDP7IPAe4B3+X/f+0fGrfDpGEWRmCYIE7p/dfVk4e5eZTct45OMb5YtQcrgUuMrMrgRGAxMIauZqzawqrI0bMZ9BNacOA+6+zt1PdPdZ7j6L4KJygbvvJHiE2c3hKNW3AgfSzQVSOczsCuBLwFXufiRj0QrgejMbZWazCQaoPFeOGCWnRmBOODqummAgyooyxyT9CPtQ/QDY6O5/k7FowEc+Snm5+xJ3nx5e664HnnL3DwBPA9eGxUbMsVNN3PD3BHAlQYf4I8AflzccyeE7wCjgl2Ft6jPu/nF3X29mDwEbCJpZP+nuXWWMU3px904z+xSwEogD/+ju68sclvTvUuAmYJ2ZrQnn/TnBIx4fMrOPEIzyv65M8UnhvgQ8YGZfA1YTJOnDnp7YICIiIhJBak4VERERiSAlcSIiIiIRpCROREREJIKUxImIiIhEkJI4ERERkQhSEiciIiISQUriRERERCJISZyIyCCZ2Q/N7D1mVmtmPzezPyx3TCIy/CmJExEZvHlAC8Gjfr7q7o+WOR4RGQH0xAYRkUEwsxhwENgL3OXu3yxzSCIyQqgmTkRkcOYAKeBDwMfNLFHecERkpFASJyIyOPOAX7r7U8CLwM1ljkdERgglcSIigzOPIHkD+AawxMyqyhiPiIwQ6hMnIiIiEkGqiRMRERGJICVxIiIiIhGkJE5EREQkgpTEiYiIiESQkjgRERGRCFISJyIiIhJBSuJEREREIkhJnIiIiEgE/X9o1Ewf1HdD8QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.signal as sig\n",
    "\n",
    "\n",
    "N = 10000  # number of samples for input signal\n",
    "K = 50  # limit for lags in ACF\n",
    "\n",
    "# generate input signal\n",
    "np.random.seed(5)  # normally distributed (zero-mean, unit-variance) white noise\n",
    "x = np.random.normal(size=N)\n",
    "# impulse response of the system\n",
    "h = np.concatenate((np.zeros(10), sig.triang(10), np.zeros(10)))\n",
    "# output signal by convolution\n",
    "y = np.convolve(h, x, mode='full')\n",
    "\n",
    "# compute correlation functions\n",
    "acfx = 1/len(x) * np.correlate(x, x, mode='full')\n",
    "acfy = 1/len(y) * np.correlate(y, y, mode='full')\n",
    "ccfyx = 1/len(y) * np.correlate(y, x, mode='full')\n",
    "\n",
    "\n",
    "def plot_correlation_function(cf):\n",
    "    cf = cf[N-K-1:N+K-1]\n",
    "    kappa = np.arange(-len(cf)//2, len(cf)//2)\n",
    "    plt.stem(kappa, cf, use_line_collection=True)\n",
    "    plt.xlabel(r'$\\kappa$')\n",
    "    plt.axis([-K, K, -0.2, 1.1*max(cf)])\n",
    "\n",
    "\n",
    "# plot ACFs and CCF\n",
    "plt.rc('figure', figsize=(10, 3))\n",
    "plt.figure()\n",
    "plot_correlation_function(acfx)\n",
    "plt.title('Estimated ACF of input signal')\n",
    "plt.ylabel(r'$\\hat{\\varphi}_{xx}[\\kappa]$')\n",
    "\n",
    "plt.figure()\n",
    "plot_correlation_function(acfy)\n",
    "plt.title('Estimated ACF of output signal')\n",
    "plt.ylabel(r'$\\hat{\\varphi}_{yy}[\\kappa]$')\n",
    "\n",
    "plt.figure()\n",
    "plot_correlation_function(ccfyx)\n",
    "plt.plot(np.arange(len(h)), h, 'g-')\n",
    "plt.title('Estimated and true impulse response')\n",
    "plt.ylabel(r'$\\hat{h}[k]$, $h[k]$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Why is the estimated CCF $\\hat{\\varphi}_{yx}[k]$ not exactly equal to the true impulse response $h[k]$ of the system?\n",
    "* What changes if you change the number of samples `N` of the input signal?\n",
    "\n",
    "Solution: The derived relations for system identification hold for the case of a wide-sense ergodic input signal of infinite duration. Since we can only numerically simulate signals of finite duration, the observed deviations are a result of the resulting statistical uncertainties. Increasing the length `N` of the input signal improves the estimate of the impulse response."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "**Copyright**\n",
    "\n",
    "This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples*."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
